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ABSTRACT 

The modelling of nonlinear and multidimensional systems from input and/or output 
measurements is considered. Tensor concepts are used to reformulate old results and develop 
several new ones. These results are verified through non-trivial computer simulations. 

A generalized tensor formulation for the modelling of discrete-time stationary nonlinear 
systems is presented. Tensor equivalents of the normal equations are derived and several efficient 
methods for their solution are discussed. Conditions are established that ensure a diagonal 
correlation tensor so that a solution can be obtained directly without matrix inversion. 

Using a tensor formulation, a new proof of the Generalized Lattice Theory is obtained. 
Tensor extensions of the Levinson and Schur algorithms are presented. 

New two-dimensional (2-D) lattice parameter models are derived. Using the tensor form of 
the Generalized Lattice Theory the 2-D multi-point error order-updates are decomposed into 
0(N 2 ) single point updates. 2-D extensions of the Levinson and Schur algorithms are given. The 
quarter plane lattice is considered in detail, first in a general form, then in forms which reduce the 
computational complexity by assuming shift-invariance. 

Based on the 2-D lattice, a new nonlinear lattice model is developed. The model is capable 
of updates in the nonlinear as well as time order. 
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I. INTRODUCTION 



Linear, time(shift)-invariant systems have been exhaustively studied and 
their properties and behaviour are well known. These systems form the 
foundation of all engineering and scientific disciplines. However, they represent 
only an approximation of reality. This fact, of course, does not diminish their 
utility. Mathematical models employing the assumptions of linearity and shift- 
invariance often provide results sufficiently accurate to be of practical use. For a 
large class of systems, however, these assumptions cannot be justified and so 
alternate models must be used. Mathematical models capable of representing 
nonlinearities and methods for their identification from system measurements are 
the major ^topics explored in this thesis. Due to the particular treatment of 
nonlinear models chosen, multidimensional linear system modelling is also 
investigated. The assumption of shift-invariance will not be relaxed. 

It will be useful to formulate a geometric framework in which to solve the 
nonlinear modelling problem. The motivation for this is simple. The transition 
from physical problems to geometric ones allows many diverse phenomena to be 
handled with a common set of mathematical tools. This relieves the burden of 
having to invent new mathematics for each new situation. Instead, the well 
understood language of geometry is used to tackle many classes of problems. 
Kron [Ref. 1: p. 197] states very clearly the rationale that allows the real 
physical problem to be converted to an equivalent geometric one in the following 
passage. 

A set of n equations with n variables (with time as a parameter) mav represent 
either the performance of a dynamical system with n degrees of freedom or the 
motion of a point along a curve located'in an n-dimensional hy pothetical space 
and expressed along some frame of reference. 

The basic approach taken in this dissertation, with respect to the modelling 
of nonlinear systems, is to represent them as a linear combination of nonlinear 
functions of the data. This allows linear algorithms to be applied in the solution 
of the modelling problem. In the process of solving the nonlinear problem several 
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new multidimensional lattice structures are developed. 

The inputs and outputs of the unknown system will be treated as random 
signals (not in general gaussian.) The problem of transforming random processes 
into geometric quantities has two solutions. These are introduced here in order 
to avoid confusion later and also in part to justify the tensor formulation that is 
used in the sequel. 

A random process X may be defined as the assignment of a function 
{x(t,w), teT}, to every outcome, w in a sample space Q. Of interest in this 
dissertation is the case when T is a discrete and finite index set. 

One way of geometrically visualizing this random process is to consider a 
Hilbert (or inner product) space, S, (in general infinite dimensional) of random 
variables, that is the vectors or elements of the space are random variables. 

^Fixing t=t 0> x(t 0 a; ) is a random variable and so is a vector in S. The random 
process, X. a time series of random variables, is a series of vectors in S. or a curve 
in S [Ref. 2: p. 27]. The components of the vectors comprising X are indexed by 
the parameter u. The required inner product on this space is defined in terms of 
the statistical correlation, ie; E{x(t 1? u;)x(t 2 ^)}. This approach has proved highly 
successful in many applications [Ref. 3]. 

An alternate approach is to consider that the random process X. consists of 
vectors in a function space. In this interpretation the random process is an 
ensemble of time functions. {x(t,u;),teT}, indexed by w. Each of these time 
functions (generally referred to as realizations) is a vector in the function space. 
There will exist a large (in general infinite) number of such vectors corresponding 
to each possible outcome, weft. The components of the vectors are indexed by 
the parameter t. There is no need to define a metric on this space. Any 
expectations that are required must be calculated over the ensemble of vectors. 

This second approach will be the one that is followed throughout this work. 
It will lead to many interesting and novel interpretations of known algorithms 
and also will be used to derive several significant new results. 

While vectors are sufficient to provide a complete characterization of 
discrete-time, one-dimensional linear systems, general nonlinear systems with 
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memory require the use of higher order geometric objects to obtain convenient 
descriptions. It is shown in this dissertation that a particular class of these 
geometric objects that extend vector concepts in a natural way and provide an 
ability to deal with nonlinearities in an organized fashion are tensors. 

The contraversial Sapir-Whorf hypothesis from linguistics [Ref. 4] states 
that the constructs of a language define the boundaries of thought. 
Mathematics is a legitimate language. It is a well defined set of rules used to 
communicate ideas. If the mathematics that is employed in the solution of a 
problem is constrained, then it is conceivable that certain solutions may not be 
arrived at, or even that the problem may remain unsolved. In Electrical 
Engineering, vector calculus and linear algebra are the major mathematical tools. 
They are adequate to explain such diverse phenomena as the propagation of 
electromagnetic waves or the behaviour of one dimensional linear systems. More 
complex problems have also been solved using this theory by forcing them to fit. 
but the notation can become awkward. Tensor analysis is a convenient 
mathematical framework in which to deal with nonlinear and multidimensional 
signal processing problems. It provides an algebra for manipulating objects of 
higher dimension than two, which is all that can easily be handled using linear 
algebra. Importantly, tensor algebra furnishes a system of notation which is 
powerful, yet compact. 

The Electrical Engineer’s experience is usually limited to ordinary. 
Euclidean geometries. Physicists, around the turn of the century, began to 
realize that other more complex types of geometries were equally valid and 
important. In fact. Einstein showed that the world we live in is neither 
Euclidean nor is it simply three dimensional. 

The arguments outlined above provide the motivation for this study of the 
utility of tensor concepts in Electrical Engineering, specifically in the area of 
discrete signal processing. Although this work examines but a fraction of the 
possible applications in this field, it proves that tensors can lead to useful results 
and that they warrant further consideration, particularly in problems involving 
spaces of higher dimensions. 
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A. HISTORICAL BACKGROUND 

Tensor analysis has evolved in this century, originating with two Italian 
mathematicians in 1900; Ricci and Levi-Civita. Many of the early contributions 
are due to Einstein who required tensor concepts in the development of his 
general theory of relativity. Recent books on the subject include Golab, Synge 
and Schild, and Young [Ref. 5,6,7], 

Kron [Ref. l] in the early 193CTs made use of tensor concepts in Electrical 
Engineering. He appears to be first to do so. His work was mainly concerned 
with the analysis and design of electrical networks and rotating machinery. Since 
that time there have been few papers that deal with tensors in the context of 
Electrical Engineering. 

Volterra [Ref. 8] laid the foundation for nonlinear system analysis in the late 
nineteenth century. He studied functionals or functions of functions. He 
proposed a series of increasing order functionals as an approximation to any other 
functional. Frechet [Ref. 9: p. 517] later showed that this series was a complete 
representation and converged uniformly. This series has since become known as 
the Volterra series. We shall study this series in detail in Chapter 3. 

The first application of the Volterra series to nonlinear systems was done by 
Norbert Wiener in the 1930’s. Wiener * also made several other significant 
contributions to nonlinear theory, such as the introduction of the Wiener G- 
functionals [Ref. 10]. They posses the property of orthogonality when the system 
input is white Gaussian noise. The two theories (Volterra and Wiener) form the 
basis of almost all significant work to date on nonlinear systems. 

One of the first practical methods of system identification was proposed by 
Lee and Schetzen [Ref. 1 1] . Their method takes advantage of the orthogonality of 
the Weiner G-functionals by employing a cross-correlation technique to identify 
system parameters. 

The study of discrete-time nonlinear systems has gained importance with the 
advent of the digital computer. It appears that the idea of a discrete Volterra 
series first appeared in the mid 1960*s (see for example [Ref. 12].) The use of 
tensor techniques in the study of nonlinear systems has received little attention. 
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Sandor and Williamson [Ref. 13] made some use of them in the study of 
continuous systems. More recently Parker and Thomas [Ref. 14] proposed the 
idea of using tensor methods for the analysis of nonlinear discrete-time systems. 
Their techniques for system identification involved the use of deterministic input 
signals to extract system parameters. 

The Volterra series is non-recursive and so a discrete form cannot represent 
an infinite memory system. This is equivalent to trying to represent an infinite 
memory linear system with a finite length impulse response. This can pose 
implementation difficulties for systems with long memories. One possible solution 
is the use of a recursive model. There has until very recently been little written 
about this because of the difficulty in analysing system stability. Parker and 
Perry [Ref. 15] have proposed a discrete nonlinear ARM A (auto-regressive 
moving-average) model, however, no stability implications were considered. Also 
Parker, Mayoral and Thomas [Ref. 16] proposed an Adaptive Kalman Identifier 
or RLS (Recursive Least Square) type algorithm for the identification of non- 
linear ARMA systems. Zarzycki and Dewilde [Ref. 17] and Zarzycki [Ref. 
18,19,20,21] have proposed a nonlinear lattice structure. Again the stability of the 
resulting models is not discussed. Some nonlinear systems are inherently 
recursive (eg: the phase locked loop) so that this remains an important area for 
research. 

Recently, several books dealing exclusively with nonlinear systems theory 
have been published. The book by Schetzen [Ref. 9] concerns itself with 
continuous systems. It provides a very thorough but readable development of the 
classical concepts. Also of interest is a short appendix outlining the history of 
nonlinear systems theory. A book by Rugh [Ref. 22]. is an important contribution 
as it includes discussions of discrete theory. 

B. DISSERTATION OVERVIEW 

Because the typical reader of this dissertation will not have a background 
which includes tensor calculus it was felt that a chapter covering some 
fundamental concepts should be included. This material was considered to be of 
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central importance to the work that followed so it remains as a chapter rather 
than being relegated to an appendix. Readers familiar with tensor concepts may 
wish to skip most of Chapter 2, although, a cursory look is recommended to 
ensure that the notation is clearly understood. 

Chapter 3 begins with a review of the traditional Volterra theory of nonlinear 
systems. Both continuous and discrete-time systems are discussed. The 
discrete-time tensor equivalent of the discrete Volterra series is deduced. An 
alternate nonlinear tensor formulation is presented along with a discussion of its 
relationship to the Volterra series. This alternate formulation will be used in 
most of the work that follows. 

Next, methods for the identification of model parameters are examined. A 
tensor equivalent of the normal or Weiner-Hopf equations is formulated. Several 
recursive in time algorithms are included as examples of the application of 
traditional linear modelling methods to the nonlinear tensor formulation. 
Nontrivial numerical simulation results are included. 

The advantages of using alternate coordinate systems are then investigated. 
It is shown that by proper choice of coordinate systems and input signals the 
identification process can be significantly simplified. The nonlinear tensor 
formulation is extended to include recursive type models. It is shown that the 
Yule-Walker equations have a tensor counterpart which can be solved for the 
model parameters. Several of the new results of this chapter have already been 
published [Ref. 23]. 

In Chapter 4 a review of modern lattice theory is presented. Although the 
results themselves are not new the approach is novel. Tensor concepts are used 
to derive the lattice filters presented by considering orthogonalizing coordinate 
transformations. Generalized forms of the Levinson and Schur algorithms are 
also presented and proven. These important algorithms are well known in linear 
matrix theory [Ref. 3] and their generalization in tensor form is a significant 
result. 

Chapter 5 breaks new ground by applying the lattice theory of Chapter 4 to 
the problem of modelling two-dimensional data fields. Simplifications due to an 
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assumption of shift-invariance are studied. Several different configurations are 
considered. Simulation results are included to prove the validity of the theory. 
Some implementational aspects of the algorithms are considered. In particular a 
systolic array is deduced for one of the two-dimensional lattice algorithms 
presented. This result demonstrates that the new algortihms are amenable to 
implementation in dedicated VLSI hardware. 

In Chapter 6 a nonlinear lattice is formulated, again based upon the theory 
presented in Chapter 4 and 5. This is a new result. The lattice structure 
proposed differs from those of previous researchers in that it is recursive, not only 
in time order, but also in nonlinear order. For example, one can obtain the 
optimal cubic model from a knowledge of the optimal quadratic model. Once 
again nontrivial simulation results are included. 

Chapter 7 is a summary of the new results presented in this dissertation. It 
draws conclusions about these results and outlines some important unanswered 
questions as possible topics of future research. 

Two appendices are included. 

Appendix A contains an alternate proof of Theorem 4.4. This proof uses the 
Hilbert space formulation described in this introduction. It is included for two 
reasons; first, to illustrate this alternate formulation and second, to provide 
additional insight into this theorem which forms the foundation of Chapters 5 
and 6. 

Appendix B contains listings of the FORTRAN programs used in the 
simulations presented in this thesis. 
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II. MATHEMATICAL BACKGROUND 



This chapter presents a brief overview of the mathematical tools that are 
used in the dissertation. It begins with a discussion of linear, bilinear and 
multilinear functionals, and it is shown that these can be represented by tensors. 
Some customary conventions which simplify notation are introduced, and some 
useful tensor operations are presented and discussed. This is meant to be an 
introduction to the subject of tensor analysis. Only those concepts that will be 
used in the remainder of the dissertation are presented. Some proofs are included 
to act as examples. Many others are not presented here, since the interested 
reader can find them in the references [Ref. 5,6,7]. The discussion assumes a 
thorough knowledge of linear algebra. 

A. LINEAR FUNCTIONALS 

We say that. V, is a vector space over a field of scalars. F. if the operations 
of scalar multiplication and vector addition are defined such that the axioms of a 
vector space are satisfied (see for example [Ref. 24,25].) 

Consider a vector space. V, over a field of scalars, F. The elements of V are 
called vectors and will be denoted by use of boldface type, viz T. If we restrict 
ourselves to spaces of finite dimension. N. then we may write a vector T as an 
N-tuple of components and denote the vector space by V N . The components of a 
vector T will be denoted by T A . where A =1....N. Writing a vector as a set of 
components implies the existence of a basis. We will denote a basis for V N as the 

set of vectors A = {a, a N ). Thus an arbitrary vector T in V N can be written 

as a linear combination of these basis vectors, viz.. 

T (2.1) 

a - 1 

In order to maintain generality it is not necessary to commit to any specific 
vector space at this point. Likewise we allow the basis to remain arbitrary. 
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The following definitions and theorems are presented essentially without 
proof. 

1. Definition 2.1: Linear Functional 



If V N is a vector space over a field of scalars F, then, a linear 



transformation, H (the reason for the boldface will become apparent shortly, (see 
eqn (2.6))), of V N into F, is known as a linear functional or linear form on V N . 
We can indicate this transformation by 



2. Theorem 2.1 

The set of all linear functionals on V N forms a vector space of the same 
dimension as V N . This space is known as the dual vector space and is indicated 
by V N . 

3. Theorem 2.2 

If V N is a vector space over the field of scalars F, with basis A = 
{a,, ..., a N }, then the set of linear functionals A = {b 1 .... b N }, (defined so that the 
A-th functional. b A . operating on an arbitrary element of V N . say T. yields the A- 
th component of T. namely T A ). form a basis for V N . The defining property can 
be expressed mathematically as 

b A (T) = T a for all A C 1 N (2.3) 



We call the functional b A the A-th coordinate function since when applied 
to a vector T it yields the A-th coordinate, namely T A . The set of these linear 



It is interesting to note that this choice of basis for the dual space leads 
to the property that 



H(T) = c where ceF and TgV n 



( 2 . 2 ) 



N 



where T - £}T A a A 



functionals. b A . A € {1 X} comprising A is known a> tin' dual basis of A. 




(2.4) 



To show this we proceed as follows: from (2.1) it follows that 
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b A (T) = b|f;ra,| (2.5a) 

= Eb A (T'aJ (2.5b) 

( 1=1 X 7 

Since b A is a linear functional (2.5b) can be written as, 

b A (T) = SrbVaJ (2.5c) 

However, from (2.3) it is known that b A (T) = T A . Thus (2.5c) implies (2.4). 

The existence of a basis for the dual space implies that any vector (linear 
functional) H in V N may be written uniquely as a linear combination of the 
elements of the dual basis. Therefore, any linear functional can be represented 
uniquely by an N-tuple of components. Thus 



H = £H A b A (2.6) 

A= 1 

From (2.6) one can write 

H(T)= £H,b A (T) (2.7a) 

Using (2.3), (2.7a) becomes 

H(T) = £H a T a (2.7b) 

A= 1 

Alternately, in matrix form, if 

H = Hj H : • • • (2.8) 

T - T 1 T 2 • • • T nt (2.9) 

then (2.6) can be written as: 

H(T) = HT. (2.10) 



We notice that the two vectors. H and T. are defined relative to different 
basis and that they belong to different vector spaces. The vector H. defined 



20 



according to the dual basis A is called a covariant vector. The vector T. defined 
according to the regular basis. A. is known as a contravariant vector. As a 
convention, whenever a subscript is used to index the components of a vector it 
is understood that the vector is being expressed according to the dual basis, A, in 
the dual vector space V N , and is a covariant vector. Similarly, when a 
superscripted index is used the components are assumed to represent a 
contravariant vector in the vector space V N according to the regular basis A. 

Equation (2.7b) represents what we normally think of as a vector inner 
product. We usually do not think of the two vectors as coming from different 
vector spaces. The reason for this is that in the familiar rectangular cartesian 
system of coordinates, the regular and dual basis are identical and so there is no 
need to differentiate between covariant and contravariant vectors. In other 
systems of coordinates the distinction must be made in order that the 
relationships have meaning. To perform a vector inner product, one vector must 
be covariant and the other must be contravariant. For example consider the 
vector T as illustrated in Figure 2.1. It has components l 3 T with respect to 
basis {e^} and components [-1 2’ T with respect to basis {e, -,e 2 '}. 

A measure of the length of vector T in the rectangular coordinate system. 
{e,.e 2 }, can be computed from the expression 



i 

I 



! T 






(2.11a) 



A= 1 




(2.11b) 



(JO)’ : (2.11c) 

The answer is correct because in the rectangular Cartesian coordinate system 
there is no need to distinguish between covariant and contravariant vectors, ie: 
T' 1 = T % . However, an expression similar to (2.11c) in the oblique coordinate 
system, {e, ,e 2 }. does not yield a measure of the length of the vector. 
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Figure 2.1: Vector T expressed according to the two bases {e,.e ; } and {e,-.e : -}. 
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(2.12a) 




= 5 1 ! (2.12b) 

The answer is incorrect since both vectors in expression (2.12) are contravariant. 
To obtain a correct answer, one of the vectors would have to be made covariant. 
This involves introduction of a metric tensor. The interested reader can find a 
discussion of this concept in [Ref. 6]. In the case of rectangular coordinate 
systems this metric tensor is the identity matrix. Our intent here was only to 
indicate that in any system other than rectangular Cartesian, strict attention 
must be paid to the character of the vectors. 

B. TENSORS 

In the previous section the concept of covariant or contravariant vectors has 
been established. In this section a more formal definition of these quantities is 
presented. 

Suppose we have N variables x 1 . x 2 . .... x N . then a set of values of these 
variables is called a point. The variables themselves are called coordinates (or 
components.) The totality of all points, as each of the variables (coordinates) 
x A , A - 1, .... N. vary over their entire specified range, constitutes an N-dirnensional 
space, denoted by V N . 

1. Definition 2.2: Contravariant Vector 

A contravariant vector T. is defined on the basis of the transformation of 
it" components upon transition from one coordinate system to another. For 
coordinate system (A) the components of T are an N-tuple of numbers designated 

HS 

T a (A 1 \) 

Cpon transition to another coordinate system (A"), if the components of T 
transform according to the rule 
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(2.13) 



T y = 



N 

E 



a= i 



8x v 

9x a 



T A 



where x A and x A define the coordinates of a point in the old (a) and new (A') 
coordinate systems, then T is said to be a contravariant vector. 

2. Definition 2.3: Covariant Vector 

A covariant vector U, is defined on the basis of the transformation of its 
components upon transition from one coordinate system to another. For 
coordinate system (a) the components of U are an N-tuple of numbers designated 
as 



U A (A = 1, .... N) 



Upon transition to another coordinate system (A'), if the components of U 
transform according to the rule 



Uv= e|^u a 

A= i d x 



(2.14) 



then U is said to be a eovariant vector. 



Q X A ' 

In equation (2.13) the quantity ~~~r represents the partial derivative of 



the new (primed) coordinates with respect to the old coordinates. Similarly, in 

equation (2.14) —jr represents the partial derivative of the old coordinates with 
d x A ' 

respect to the new. primed, coordinates. In general these quantities can be 
arranged into a two-dimensional matrix of numbers. 

3. Example 2.1 

Consider the following parametric description of a curve: 



d - f,(t) 



(2.15a) 



x 2 =fj(t) (2.15b) 

x* - f*(t) (2.15c) 

We consider the three quantities x 1 . x 2 . x 3 , to be the coordinates of a point or 
equivalently the components of a vector in a three dimensional space. We leave 
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the basis arbitrary. Indeed the equations we write will be true regardless of the 
choice of basis. This is the inherent advantage of tensor analysis since it allows 
expressions to be written which are invariant with respect to the coordinate 
system. 

We can now define new, primed, coordinates as functions of the old 
coordinates. For example, if we arbitrarily choose the following coordinate 
transformation 



x 



r 



x 



2 ' 




(2.16a) 



(2.16b) 



x r = g 3 (x\x 2 .x 3 ) 




\ 



then, the quantities 



3x a ' 
9x a 



can be written in matrix form as 





(2.16c) 



(2.17) 



According to equation 2.13 any contravariant vector, say T. with components T ; ‘. 
can be expressed in the new (primed) coordinate system as 

ryX ' _ o x r-pA 

^ Rv A 

X~ 1 



or in terms of components 




(2.1Sa) 



(2.18b) 
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(2.18c) 



rp3 ' 



— T 1 



OT 2 





T 1 ) 



4. Definition 2.4: Contravariant Tensor of Order 2 

A set of N 2 numbers T A ", where A and // = are said to be the 

components of a second order contravariant tensor if, upon transition to another 
coordinate system, they transform according to the rule 



N N A ’ r\ * 

x=\ „=i 9xA 9x " 



(2.19) 



5. Definition 2.5: Covariant Tensor of Order 2 



A set of N 2 numbers lfi„. where A and /a = 1,....N are said to be the 
components of a second order covariant tensor if, upon transition to another 
coordinate system, they transform according to the rule 

\ 

( 2 . 20 ) 



N N 



u 



a>-= E E 



9x-* dx" 



A=1 M =1 



8x A dx M 



Similarly, tensors of higher order can also be defined. In the general 
case, it will no longer be possible to use different letters to denote indices. In this 
case indices with sub-indices will be utilized, namely A 1} A 2 , .... A N . 

6. Definition 2.6: Contravariant Tensor of Order p 

A set of N p numbers T* 1 \ where X i = 1, N for i = l,...,p. are said to be 

the components of a p-th order contravariant tensor if. upon transition to 
another coordinate system, they transform according to the rule 



T V 



N 

V 

A, j 



E 



a. 9 



A, ' 



d x 



8x V 

9x Ap 



( 2 . 21 ) 



7. Definition 2.7: Covariant Tensor of Order p 

A set of N p numbers 1 x . where A 1 N for i l. ,p. are said to be 

the components of a p-th order covariant tensor if. upon transition to another 
coordinate svstern. thev transform according to the rule 



Ua. 



V 



E • ■ E Ua, 



A . = 1 



V 



9.y A ‘ 

Ap 9x A * 



_9aA 

9x Ap 



( 2 . 22 ) 
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8. Definition 2.8: Mixed Tensor of Order p 

A set of N p numbers S'* 1 \ > where A. = 1 N for i = 1 p. are said 

q + 1 p 

to be the components of a p-th order mixed tensor, with q contravariant and (p- 
q) covariant indices if, upon transition to another coordinate system, they 
transform according to the rule 






9x A 1 



9 x A q 9 x Aq+1 
9 x A ’ 9 x A q+1 




(2.23) 



We have already seen two examples of mixed tensors. The first is the 
Kronecker delta 



6\ = 



1 whenA - n 



(2.24) 



0 whenA # (i 

To see that the Kronecker delta is in fact a mixed tensor we must test to see if it 
transforms according to the rule given in equation (2.23). We must prove that 
the following relation is true 



r A ' ip yi 9x J 9^ p 

* ~ E E a A 8 * 

A = 1 U = 1 9X 9X 



(2.25) 



We begin with the right hand side of (2.25) 



V V — — 

2-J 2-J r A CX U 

a = i „ = i 9 x 9x" 



E 



N 9x A ' A ,a Sx" 



A- 1 9x A 



E s ' 



dx p 



1 2.26a) 



y, 9 x A 9 x a 
9 x a 9 x" 



(2.2Gb) 



dx 

dx'' 



A' 



!.26c 



- (2.26d) 

Therefore, we conclude that relation (2.25) holds and so the Kronecker delta is in 
fact a second order tensor of mixed character. 
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The other mixed character tensor that we have already worked with is 
the one appearing in the formulae of transformation (definitions 2.2 through 2.8). 
that is the the partial derivative of the new coordinates with respect to the old 
(and the old with respect to the new). We will not prove that this is in fact a 
tensor of the type stated although the proof is straight forward. It is instructive 

3 x^ 3 x*^ • 

to note, however, that the two quantities — — and r are inverses of each 

3x A 3x 

other. 

As a final note, vectors are tensors of order one. Also scalars are 
considered to be tensors of order zero. They are sometimes called invariants 
since their representation is independent of the coordinate systems used. 

C. NOTATIONAL CONVENTIONS 

There are two widely accepted conventions that simplify notation and 
unquestionably save much writing. The first is known as the summation 
convention. Historically, it was first used by Einstein. He noticed that in almost 
all cases there is really no need to explicitly write summation symbols. 
Summation can be implied whenever an index is repeated in an expression, once 
as a superscript and once as a subscript. The repeated index is allowed to take 
on all permissible values and the resulting terms are summed together. This type 
of index is often referred to as a dummy index. But what are permissible values 
for the index? This question leads to the second convention. Normally, the range 
of the index will not be explicitly stated. By convention it is understood that all 
greek subscripts and superscripts appearing in an expression will take on all 
values from 1 to N. where N is the dimension of the vector space in which we are 
working. In later chapters we will find it more convenient to allow indices to run 
from 0 to N. The dimensionality of the vector space will thus be N-t-1. An 
additional convention which we shall find useful is to reserve latin indices to 
indicate that we are dealing with a particular component of a quantity. In most 
books this is indicated by surrounding the particular index with parenthesis. 
However, we will reserve parenthesis to indicate exponentiation. The conventions 
adopted here will be used throughout the sequel. In exceptional cases, where 
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some deviation from them is required, we will explicitly state the meaning of the 
notation. As an example of the use of these conventions consider the expression 

Y a = E H a „T" for A = 1. ..., N (2.27) 

It can be written more succinctly as: 

Y a = H a „T" (2.28) 

Another convention which has already been used is now formally introduced. 
Every tensor quantity will be given a distinct base letter. Upon a change of 
coordinates, the base letter will be maintained in order to indicate that the 
quantity has not been modified, only the representation has changed. The 
coordinate system used is indicated by the sub- or superscript used to index the 
components of the quantity. We therefore, will refer to different coordinate 
systems simply by the index letter used to indicate the components. For 
example, a vector T has components T A in the (a) set of coordinates, while it has 
components T A ' in the (A') coordinate system. We note that using this 
representation, scalars appear identical in all coordinate systems, which is 
desirable. 

1. Example 2.2 

The following example serves not only as an illustration of the 
conventions presented in this section, but also as a concrete (and presumably a 
somewhat familiar) illustration of the two types of vector. Consider an invariant 
function of the coordinates, f f(xh x 2 . x 3 ). The differential of this function is given 
by 



df 



0 f 
dx 1 



dx’ 



df 

dx J 



dx 2 




(2.29) 



We can consider dx A to be the components of a contravariant vector representing 
an infinitesimal displacement expressed according to some baW A = {a,. a ; a,}. 
We can write this as: 



29 



(2.30) 



dx = 2 dx>a A 

A= 1 



We have already shown that components of a contravariant vector transform 
according to 

dx A ' = — — sr-dx A (2.31) 

9x a v ’ 

The gradient is also a vector whose components are given by: 

vr, = (2.32) 



Upon a change of coordinates the values of the new components, can be 

deduced by application of the chain rule; 

9f 9 x a 



vU- = 



9x A 9x A 



(2.33) 



It is clear that the components of the gradient transform according to the rule 
given in equation (2.14). The gradient must, therefore, be considered to be a 
covariant vector. 

2. Example 2.3 

Although it has already been stated (section 2.1) that linear functionals 
can be considered to be covariant vectors, this fact has not been proven. In this 
example we will show that any linear functional, say H. with components H A that 
transforms an arbitrary contravariant vector, say T. (according to equation 
(2.7b)) to yield an invariant, satisfies the definition of a covariant vector 
(equation (2.14)). Equation (2.7b). which defines a vector inner product is 
repeated here for convenience. 



H(T) - IfiT A 



(2.34a) 



The quantities T A are the components of an arbitrary contravariant vector. 
Because of the assumed invariance we may write 



H a T a = H a -T a 

From the definition of a contravariant vector, (equation (2.13)) 



(2.34b) 
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(2.34c) 




Equation (2.34b) becomes 




(2.34d) 



Rearranging yields, 




(2.34e) 



Equation (2.34e) implies 




(2.34f) 



since the contravariant vector, T, was arbitrary. A simple change of variables in 
this last expression yields 



thus justified in calling linear functionals covariant vectors. 

D. BILINEAR AND MULTILINEAR FORMS 

We next consider higher order functionals. In particular we will start with 
the bilinear form or bilinear functional. 

1. Definition 2.9: Bilinear Functional 

A mapping f of a pair of vectors, say T e U N and S e V M into a field of 
scalars, is a bilinear functional or bilinear form if f(T,S) is a linear function of 
T and S taken independently. We will only consider cases when N = M and 
U N V N . 

Once again choosing an arbitrary basis A = {a,. ... a N } we may express 
two arbitrary contravariant vectors as linear combinations of these basis vectors. 

T = T A a^, and S - S^a^ 

The bilinear functional f can then be written 




(2.34g) 



which is identical to the equation defining a covariant vector. (2.14). We are 
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f(S.T) - f(S*a M; T A a A ) 



(2.35a) 



= S*f(a^T A a A ) (2.35b) 

= S M T A f(a M ,a A ) (2.35c) 

The bilinear form is thus completely determined by the N 2 quantities f(a^,a A ). We 
will write these components of f as f^ A . Using this shorthand, equation (2.35c) can 
be written as 

f(S,T) = f; A S*T A (2.35d) 

In matrix notation the bilinear form takes on the familiar appearance 

f(S,T) = S t FT (2.36) 

where F = [f |lA ] 

If the two vectors S and T are equal then the bilinear form reduces to the 
well known quadratic form 

f(T,T) = f MA T*T A (2.37a) 

or in matrix notation 

f(T.T) T t FT (2.37b) 

We will be interested in the behaviour of the components. f^ A , of the 
bilinear form. f. upon transition from one coordinate system to another. We 
establish their tensor character in the following example. 

2. Example 2.4 

In Example 2.3 we showed that a linear functional satisfied t lie definition 
of a covariant vector. Here we will show that a bilinear functional satisfies the 
definition of a covariant tensor of second order. It h necessary for the discussion 
that follows in later chapters to establish the tensor character of bilinear 
functionals. Since the bilinear functional yields an invariant (scalar), we may 
write 

= f^-S" T v (2.38a) 



32 



(2.38b) 



= WS" 



9x*' 

9x p 




dx M ' 

9x" 



9x A 

9 x a 



Va-S"T a 



Rearranging yields 



(U> 



9x M 

9x^ 



^Lf^S'T* = 0 

OX 



(2.38c) 



(2.38d) 



Since the contravariant vectors S and T are arbitrary the quantity inside the 
parenthesis must vanish. This yields the relation 



dx* 9x A 
3x"' 3 x a 



(2.38e) 



This last equation is identical to the definition of a second order covariant tensor 
(equation (2.20).) We have thus proven that bilinear functionals are covariant 
tensors of order 2. 

In general we can have m-linear functionals which map m contravariant 
vectors. T(l), .... T(m), into a scalar and are linear functions of each of the m input 
vectors taken separately. They can be written as 

f(T(l) T(m)) = T A, ( 1 ) • T A -(m)f Ai ( 2 . 39 ) 

Using identical arguments to the ones presented for linear and bilinear 
functionals, we can show that multilinear functionals are also covariant tensors. 



E. TENSOR OPERATIONS 

There are a few tensor operations that will be of considerable importance in 
later discussions. Although some have already been used we will formally define 
them here before proceeding. Only those operations that will be used in the 
sequel are presented. Others are possible and are discussed at length in the 
references (see for example [Refs. 1.5. 6. 7].) 
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1. Definition 2.10: Tensor Outer Product 



Given the components of two tensors 



r p 1 n 

r u 



and S A " +1 Ap 



^ m+ 1 ^ q 



the N* p+q * numbers R A ‘ Xp lli . . . „ given by 



_ A. • • • A„ 

R * 



_ 'T'^l c^n+I 

1 *1 ^ m+ 1 • 



(2.40) 



(2.41) 



are components of a tensor of order p + q. The operation implied in the above is 
known as the tensor outer product or simply the tensor product. 

2. Example 2.5 

As an example of the tensor product operation consider two vectors 
defined as: 



T = |T A ' and S = [S M ] (2.42) 

The tensor product is given by (equation (2.41)) 

r a " = T A S" (2.43) 

In this case the components of the tensor product can be arranged as a matrix of 
N= numbers. For simplicity consider the case when N=3. In matrix notation 

R = TS t (2.44a) 



or 



R a " = T a S" 



(2.44b) 



T’S 1 


T‘S 2 


T‘S 3 


T 2 S‘ 


T“S* 


T 2 S? . 


T 3 S' 


t ! s : 


t 3 s 3 



(2.44c) 



3. Definition 2.11: Contraction 

The operation of .netting two indices, one lower and the other upper, 
equal and summing the result is known as contraction. The result is a tensor 
which has the character indicated by the remaining indices. The contraction of a 
tensor of order pH-1 over two indices results in a tensor of order p - 1. 
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4. Example 2.6 

We may contract the tensor 



n x 



H 1 , h 1 , 

H 2 u2 2 
l n 2 n j 

H S it 3 tt 3 
1 n 2 ^ 3 



over the indices A and resulting in 



H a a = H\ + H 2 2 + H S j 



(2.45) 



(2.46a) 



= Trace|H A J (2.46b) 

a scalar which represents the trace of the matrix. 

We can consider a higher order example. Assuming a three dimensional 
vector space, consider contracting a tensor U ^ 7 over the indices A and n . The 
result will be a vector wdiose components, U 7 , are given by 

U 1 = U 1 , 1 -s UV + U *, 1 (2.47a) 

U a = U 1 , 2 + UV + U*s 2 (2.47b) 



U* = U 



1 9 
J 

] 



uV - uV 



(2.47c) 



5. Definition 2.12: Tensor Inner Product 

Suppose that after we form the outer product of two arbitrary tensors T 

and S. with components T A ' \ „ and S^"* 1 \ „ . we set the indices A 

M J ^ m ^ m A 1 ^ q 

and equal. This implies a contraction operation. The outer product h 
written (according to equation (2.41)) as 



r A| 











u 



" q 



(2.48a) 



The contraction operation yields. ( definition 2.11) 



R A| 



4 J I** j- 1 



R A ' 



“j i ' j- i 



(2.4Sb) 



It can be shown that the object. R. given in the last equation is a tensor of the 
character shown. If the original tensors. T and S. were of order m+n and p+q- 
(m+n) respectively, then the resulting tensor. R. will be of order (p-fq-2). The 
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total operation, consisting of contracting the result of a tensor product over one 
(or more) pairs of indices, each member of the pair having come from a different 
tensor, is known as a tensor inner product, or simply inner product. This 
operation is also sometimes referred to as transvection (see for example [Ref. 

5 ]) 

This operation has been used previously in the discussion of linear, 
bilinear and multilinear functionals, (see equations (2.7b), (2.35d) and (2.39) 
respectively.) It will be particularly valuable in future discussions. 

6. Example 2.7 

Some insight into the inner product operation may be gained by 
explicitly performing the two steps described above for the simple case of the 
linear functional (equation (2.7b).) 

y = H a T^ (2.49) 



We can first form the outer product by replacing one of the A indices appearing 
on the right hand side with a different letter, say p.. We may then write: 





H,T' 


H,T 2 


H,T 3 


H A T" : - 


H;T‘ 


IhT 2 


H,T 3 




H.T 1 


H.T 2 


H.T 8 



(2.50) 



The result is now contracted by equating A and p and performing the implied 
summation. 



y = H A T A 



(2.51a) 



= H,T‘ - H 2 T- - H,T‘ (2.51b) 

It is often useful to perforin both these steps (tensor outer product and 
contraction) when calculating an inner product, particularly when dealing with 
higher order tensor^. It may otherwise be difficult to keep track of how terms are 
combined or even which terms will be present in the final expression. 
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III. NONLINEAR SYSTEMS THEORY 



In this chapter several models of nonlinear systems are described. The 
discussion begins with an overview of the classical continuous Volterra series. 
Several interesting aspects of the series are examined. Next, a discrete-time 
version of the Volterra series is introduced and is used to develop an equivalent 
tensor form. Finally, a new discrete-time, nonlinear tensor model is presented 
and its relationship to the Volterra series is discussed. 

A. CONTINUOUS NONLINEAR SYSTEMS MODELS 

Traditionally, non-linear systems have been modelled using the continuous 
Volterra series expansion [Ref. 26: p. 1559]. In its most general form this can be 
written as: 

>'(') = holM 




).x(r j )dr i 



DC OC 

- f f h 2 (t,r ,,r 2 )x(r ,)x(r 2 )dr,dr 2 

— 00— OC 



(3.1) 



where x(t) is the system input and y(t) is the system output. The parameter- 
h 0 (t). hj(t.7-|). h 2 (t.r j.r 2 ). are known as the Volterra Kernels. As we will only 

be concerned with time-invariant systems, the kernels will only be functions of 
the time difference, (t r). and not the actual time. The kernels then become: 
h 0 , h, (t — j- j). h 2 (t r ,,t - r 2 ). Simple changes of variables then allow the series. 

(3.1). to be written in the form 
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y(t) = 4 



oc 




— oo 



oo oo 




+ * * * 



(3.2) 



We note several things about the expansion. First, an infinite number of 



h 0 . hj(r ,) h 2 (r 1 ,r 2 ), ••• correspond to the constant, linear, quadratic, ... terms of 
the expansion, respectively. The familiar linear system appears as a special case 
of this more general expansion (ie; the case when all kernels except h are 
zero.) The third term: 



is a bilinear term, that is it is linear in each variable x(t-rj) and x(t-r 2 ) taken 
independently (ie: assume the other variable is a constant.) In general, the 
(i+l)-st term is i-linear. It is linear in each of the i variables 
x(i -r x(t-r 2 ), x(t - rj taken one at a time. Lastly, notice that the Yolterra 
expansion is not orthogonal in the sense that the identification of the ii-th order 
kernel depends on the values of all the other kernels. They cannot be identified 
independently. 

The non-linear model represented by equation (3.1) can be visualized as 
diown in Figure 3.1. As can be seen it corresponds to a parallel connection of 
subsystems of increasing non-linearity. Each of the subsystems is homogeneous 
(except for the zeroth order subsystem) in the senile that increasing (multiplying) 
the input by a factor k results in the output of the p-th subsystem being 
increased by a factor k p . This can easily be understood by examining the 
expression for the p-th term. 



terms are required to represent the most general case. Second, the kernels 



oc oc 




(3.3) 
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(3.4) 



OC 00 



f • • J h p (r 1 , r p )x(t-7-,) ••• x(t-J- p )dr, • dr p 

-00 — 00 

Replacing x(t) by kx(t) yields 

oo oo 

f • • • f hp( r i» •••» r p) kx ( t - r i) • • • kx(t-r p )dr, • ■ dr, 



-00 —oo 



(3.5a) 



= kP f ' ' ' f b p( r i> ■■■’ r p) x ( t ~ r i) x ( t - r P ) dr i dr p (3.5b) 

— oo —oo 

The presence of a constant term in the Volterra expansion should not be 
unexpected. Consider, for example, a system whose response is given by. 

y(t) = x(t)-h (3.6) 

It is easily shown that this system does not obey the principle of superposition 
and so cannot be considered linear. This necessitates the inclusion of a constant 
term in the Volterra expansion in order to handle the general case. The usual 
procedure adopted in linear analysis, if a constant term appears, is to define a 
new output function which is the actual output less the constant term. This new 
output function is then identified in the usual fashion. The constant term must 
be separately identified. If we admit that the system is non-linear then this will 
no longer be necessary. 

There are many other aspects of continuous-time nonlinear system modelling 
that have not been discussed. In the next section discrete-time nonlinear systems 
are introduced. Many of the comments that will be made there are equally 
applicable to continuous-time nonlinear systems. However, we make them in the 
context of discrete-time, since that will make them more applicable to the sequel. 

B. DISC'RETE-TIME NONLINEAR SYSTEM MODELS 

A discrete version of the time-invariant Volterra expansion can be written as: 

y(k) - h 0 

+ E h,(A,)x(k-A,) 

q=o 
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Figure 3.1: Nonlinear Volterra Series Model 
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(3.7) 



“ S X h 2 (-\i-Ai)x(k - AJx(k A 2 ) 
a, = o a 2 = 0 



where we have assumed that the system is causal. 

For a large class of systems the summations can be truncated to N+l terms. 
Certainly truncation is required for computer implementation. This truncation 
implies that only finite memory nonlinear systems can be represented. Equation 
(3.7) is an extension of the linear Moving Average (MA) type model. In fact the 
linear model is a special case of equation (3.7). This expansion is non-recursive in 
nature, that is. it expands the present output only in terms of the present and 
past input. Past outputs are not used. 

The representation of the kernels in both the continuous and discrete forms 
of the Volterra expansion is not unique. There are. however, several special forms 
which are important. Consider a second order kernel for which 
h 2 (A,,A 2 ) = h 2 (A 2 .A,). This kernel is symmetric with respect to the two parameters A, 
and A 2 . It turns out that the kernel can always be symmetrized with no loss of 
generality. For the p-th order kernel the procedure for obtaining the symmetric 
kernel from an asymmetric one is given by [Ref. 22]: 

h,ymUi ^ 2 ) - ~rX]hp(A*(i). •••> ^»(p)) (3-8) 

n ’ r() 

where the summation is over all n! possible permutations of the p A 's. Although 
the symmetric form may contain more terms than an asymmetric form it is of 
importance because it is unique [Ref. 9: p. 43]. There can be many equivalent 
asymmetric forms of the kernel which all lead to the same symmetric kernel 
(through equation (3.S)l. The symmetric kernel thus provides a standard form 
which can be used as a reference. 

Other forms of interest are also possible. The symmetry of the symmetric 
kernel implies redundancy. This redundancy can be eliminated by use of a 
triangular kernel. Consider a kernel defined so that h tr ,( A , A p ) = 0 whenever 
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A S >A, for s < t. The domain of a second order triangular kernel is illustrated in 
Figure 3.2a. For comparison, the equivalent symmetric kernel’s domain is 
illustrated in Figure 3.2b. 

Using the triangular kernel the output of the p-th order nonlinear subsystem 
is given by 

y(k) = EE''' Eh t ri(Ai, ..., A p )x(k- A,) ■ ■ x(k-A p ) (3.10) 

V° A 2=° V° 

Notice that the limits of the summations reflect the triangular domain. This 
implies that fewer terms are included in the summations resulting in 
computational savings. The triangular kernel defined above, and used in 
equation (3.10). is not unique. Other triangular kernels can be formed by 
choosing alternate triangular domains. For example, the domain illustrated in 
Figure 3.3 could equivalently have been used. This choice corresponds to setting 
h tri ( A i, .... A p ) - 0 whenever A s > A t . for s > t. 

The output of nonrecursive models of the type presented in equation (3.7) is 
stable if the input is bounded and if the series is truncated to a finite number of 
summations. Consider an input x(n) which is bounded to be less than some 
constant M. If the series is truncated to p + 1 terms then, in the worst case the 
output will be 

y(n) $1 h -t- M E MA.) f M 2 E E MA,.A 2 ) - 

a,= o a, = o A : =o 

-m p E E h p( A i A P ) (3.n) 

A.-O -p'O 

So. a- long as the summations are bounded (which will generally be the case), the 
output will always remain finite. This guaranteed stability makes MA type 
models very attractive. As mentioned earlier, their shortcoming is their inability 
to accurately model infinite memory systems without using a large number of 
terms. 
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a. Second Order Triangular Kernel Domain 




b. Second Order Symmetric Kernel Domain 



Figure 3.2: Triangular and Symmetric Yolterra Kernel Domains 
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1. Tensor Formulation 



In order to adopt a tensor formulation of the problem we notice that 
equation (3.7) can be considered to consist of a series of increasing order 
functionals. As has been shown, these functionals can be expressed as tensors. 
We can therefore rewrite equation (3.7) as a tensor equation. 

y(k) = H 



H 



A X 



1 



Hi 



A x 
ri 



+ * * * (3.12) 

where we have defined the contravariant input vector as 

x = ,x(k) x(k-l) * • * x(k — A) • • • x(k— N )] t (3.13a) 

= x AlT (3.13b) 

This choice for x has the effect of truncating the series to N+l terms. The 
symbols H, H Ai> * • • represent the components of covariant . tensors of order 

0. 1, 2, etc., respectively. 

Examining equation (3.12) we make particular note of the following two 

aspects: 

(1) The dimension of the vector space is related to the memory order of the 
system. The vector x has N + l components implying that the nonlinear 
system contains no more than N delays. This fact is explicit in the way the 
input vectors have been defined. 

(2) The nonlinearity of the expansion is provided implicitly by the tensor outer 
product operation. It is the outer product of the vector x with itself that 
makes the higher order (2 and up) terms nonlinear. 
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Figure 3.3: Alternate Second Order Triangular Volterra Kernel Domain 
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These observations provoke speculation about the possibility of an 
alternate formulation where the roles played by the dimensionality of the vector 
space and the outer product operation are interchanged. 

2. Alternate Tensor Formulation 

Consider a parametric description of a curve in V p+1 , a p+1 dimensional 
vector space, given by; 



x°(k) = [x(k)|l°> 
x*(k) = [x(k)]W 



x(k) = 



x p (k) = [x(k)]< p ) 



(3.14) 



where the superscript in parenthesis indicates exponentiation. Given any x(k). 
the components x A (k), A = 0,...,p define a point in V p+1 . As x(k) varies with k. the 
vector x(k) will describe a curve in V p+1 . We define the components of another 
vector in V p+1 as 



x A, (k-l) = [x(k-l)j (A,) (3.15) 

Similarly, the N-th such vector can be defined as. 

x AN (k-N) = ;x(k-N)! (AN) (3.16) 

The vectors in equations (3.14). (3.15). and (3.16) will be referred to as 
observation vectors. Although, at this point they only depend on past and 
present system input values, later, in Section D, they will be generalized to 
include past outputs as well. The input and output measurements represent the 
system observations, hence the name. 

Using the theory developed in Chapter 2. concerning multilinear 
functionals, the following mathematical relation is proposed as a model of a finite 
order, finite memory nonlinear system: 

y(k) = x A °(k) • ■ • x AN (k-N)H v An (3.17) 
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where H Ao . A is an (p+l)-order covariant tensor representing a (p+l)-linear 
functional. This covariant tensor performs a role similar to that of a Volterra 
kernel. 

This model is equivalent to the p+1 term Volterra series (equation (3.7) 
or (3.12)), although it does contain many additional terms. It has the advantage 
of notational simplicity. It replaces a p + 1 term series with a single term and 
requires the specification of only a single composite kernel, H A . . . A . As will be 
shown in Section D, equation (3.17) is considerably more general than equation 
(3.12). It will be shown that Wiener type models can be obtained from equation 
(3.17) by a simple coordinate transformation. Other choices of observation 
vectors yield Autoregressive (AR) or even Autoregressive-Moving Average 
(ARMA) type models. 

In order to illustrate the correspondence of equation (3.17) to the 
standard Volterra type series expansion of equation (3.12) we present a simple 
example. 

3. Example 3.1 

Consider the truncated Volterra series expansion corresponding to 
equation (3.12): 

y(k) = H + H,/ 1 + H v xV 2 (3.18) 



where A, = 0.1 and A 2 = 0,1. and 





x 1 




x(k) 


X — 


x 2 




x(k-l) 



We may explicitly write out all the terms implied in equation (3.18). This yields 



y(k) = H 



- H 0 x(k) -r H,x(k 1) 



+ H 00 x(k)x(k) -r H 01 x(k)x(k-1) 



+ H 10 x(k— l)x(k) + H u x(k— l)x(k-l) 



(3.19) 



Now consider a model of similar order given by equation (3.17). 



y(k) - H v x A °(k)x A *(k-l) (3.20) 

where A 0 , A, = 0,1,2 and 

x>) = |x(k)] (Ao) 



x A ‘(k-l) = [x(k-l)] (A,) 

The tensor product, x A °(k)x A ‘(k-l), appearing in equation (3.20) results in a 
contravariant tensor of second order. Its elements are 



|x A °(k)x Al (k— I)) = 



1 x(k-l) 

x(k) x(k)x(k- 1) 

x< 2 )(k) x (2) (k)x(k-l 



x (2) (k- 1) 
x(k)x< 2) (k-l) 
x< 2 )(k)x< 2 >(k-l) 



(3.21) 



We notice that all the elements on or above the southwest- northeast diagonal 
are included in equation (3.19). the terms below this diagonal are not. This new 
form has not discarded any terms present in the latter version and so cannot 
represent a loss of generality. The extra terms that are included in this new form 
do not pose any significant problem. If the system does not contain terms 
involving these particular elements then the corresponding term in the kernel will 
go to zero during the identification process. Certainly in some classes of systems 
these additional elements will be important. 

C. DISCRETE NONLINEAR SYSTEM IDENTIFICATION 

In Section B a tensor formulation of the discrete nonlinear modelling problem 
was introduced. In this section methods for kernel identification are presented. 
All the methods that are discussed are statistical in nature and utilize a least- 
squares approach of parameter estimation. It will be shown that familiar 
methods used in linear systems modelling can be extended to handle the 
nonlinear case. In the first section a tensor equivalent of the normal equations, 
which can be solved for the unknown system parameters, is derived. Several 
recursive (in time) solutions to the problem are then presented. Finally, 
simulation results are offered to illustrate the effectiveness of the algorithms. 
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1. Derivation of Normal Equations 

In Section B (equation (3.17)). the following tensor model for nonlinear 
systems was proposed: 

y(k) = x A °(k) • • • x AN (k-N)H v An (3.22) 

where y(k) represents an estimate of the system output and the x A ‘(k-i)'s are 
components of contravariant observation vectors, defined in equations (3.14) 
through (3.16). 

Consider the analysis model of Figure 3.4. This diagram represents a 
conceptualization of the system identification process. The assumption is made 
that the unknown system can be represented by an equation identical to the 
model equation, (3.22), where the system parameters H Ao . ,. A are unknown. The 

parameters of the model are adjusted to best match the actual system 
parameters. A convenient measure of how well the model represents the actual 
system is the mean-square error or MSE. The MSE is a quadratic function of the 
model parameters which implies that there exists a unique minimum, or optimal 
solution, to the problem. Minimizing the MSE yields a linear set of simultaneous 
equations which can be solved for the unknown model parameters. In addition, 
the quadratic nature of the MSE allows steepest descent type, adaptive 
algorithms to be used. 

The error signal, e(k). defined as the difference between the actual 
nonlinear system output, y(k). and the output of the model. y(k), is given by 



e(k)= y(k)-y(k) 

The mean-square value of this error signal is given by 

E|e-(k)j = i:| V ( k ) y(k) 4 



(3.23) 



(3.24a) 




y(k)-x A| '(k) 



x N (k-N)H; 




(3.24b) 
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Figure 3.4: System Analysis Model 
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where E{ } indicates expectation. The optimal set of parameters can be found by 
minimizing the MSE, E{e 2 (k)}. with respect to the parameters H A( A . The 

gradient of the MSE is formed by differentiating with respect to the unknown 
parameters. 



9Ee 2 (k) 



9H V h, 



= 2Ej[y(k) - x» ' • ' x^(k-N)H Ao . [-x'‘°(k) • • • x" N (k-N)] 



Setting this equal to zero yields 



2E<’[y(k)x Mo (k)...x" N (k-N) - x~°(k)...x^(k — N)x*°(k)...x* N (k — N)H V Afj 






(3.25a) 



= 0 (3.25b) 

or. 

E{x''°(k)...x AN (k-N)x' I °(k)...x < ' N (k-N)[H Ao An = Ejy(k)x"°(k)...x'‘ N (k-N) 

(3.25c) 

The expression. (3.25c). is a tensor equivalent of the Wiener-Hopf or 
normal equations. As will be shown in Section D. these equations can also be 
used to represent the Yule-Walker equations if a different choice is made for the 
observation vectors, x(k-i). The term on the left-hand-side of equation (3.25c) is 
a nonlinear extension of the autocorrelation matrix. This contravariant tensor of 
order 2 (X + 1) includes various higher order correlations as well as the second 
order ones which arise in the linear case. The expectation on the right-hand-side 
represents a cross correlation between the output and various linear and 
nonlinear functions of the input. 

The minimum MSE can be determined by substituting (3.25c) into 
(3.24b). This yields 
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2E y(k)x °(kj • • • x N (k X)H ; 




Hi 



H„ 




A °(k)...x AN (k-N)x"°(k).. 




(3.26a) 



= E|y 2 (k)J - E|y(k)x A °(k) • • • x A »(k- N)|h a# .. (3.26b) 

Equation (3.25c) represents (p-sl) (N+1) equations in as many unknowns, 
and so can theoretically be solved for the unknown parameters H Aq . .. A . 
However, in practice the number of computations required to perform the matrix 
inversion becomes unwieldly. An nxn matrix inversion takes on the order of 0(n 3 ) 
operations [Ref. 27: p. 58], therefore, 0( p+ij 3(N ' Tl )) operations will be required in 
the solution of (3.25c). In order to make the task manageable, alternate 
algorithms that avoid direct matrix inversion, must be employed to solve the 
normal equations. 

2. The Least Mean Square (LMS) Algorithm 

The Least Mean Square (LMS) algorithm has successfully been employed 
in the solution of a variety of linear problems Ref. 28]. It is a recursive 
algorithm based on a gradient, steepest descent type of strategy. The mean- 
square error is a hyperparaboloid which is concave upward. Steepest descent 
algorithms strive to descend down this quadratic surface, towards the minimum, 
by making adjustments to the unknown parameters proportional to the gradient. 
This can be expressed mathematically as 

H Ao A N (k-l) - H Aq AN (k) - „v v -a n 00 (3.27) 

where A Jk) is the value of the model parameters at the k-th time instant 

and // is a parameter which controls the convergence of the algorithm. The 
symbol. * (k). is used to represent the gradient. 

The actual value of the gradient is given, according to equation (3.25a). 
by 
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(3.28) 



Va 0 j N <k) = 2E|e(k) x A °(k) • x^k-N) 

The LMS algorithm approximates the MSE by the instantaneous value of the 
error squared. The approximate value of the gradient is 

VA 0 ...A N (k) = - 2e(k)[x A °(k) • • • x >N (k-N)] (3.29) 

where the circumflex indicates an estimated value. Using this approximation in 
equation (3.25) yields 

Ha o , N (k-l)=H v , N (k) + 2^e(k)x A °(k) • ' • x" N (k-N) (3.30) 

Equation (3.30) gives a straight forward method of determining the 
system parameters. It involves no matrix inversion nor does it require the 
correlation tensor to be known. These two properties make the calculations 
required at each iteration very simple, so that it is possible to perform them in 
real time. It can be shown that the LMS algorithm converges to the optimal 
solution [Ref. 29: p. 578], 

For the linear case, the algorithm will converge as long as the parameter. 
fi. is chosen to satisfy. r Ref. 29: p. 578’. 

°< v < -T- — (3.31) 

max 

where A max is the largest eigenvalue of the correlation matrix. Since the 
correlation matrix is positive definite, all its eigenvalues are greater than zero. 
Therefore, the trace of the correlation matrix will always be greater than the 
largest eigenvalue. The following condition will, therefore, ensure stability. 



0 < /j < 

" Tr correlation matrix 

For the nonlinear case this translates to 



0 «-‘ fj - 



£ i: E { x A "(k)...x^(k X)x A "(k)...x K (k-X) } 
0 A 0 V / 



/ O o O \ 
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In practice a value considerably smaller than the maximum permitted by 
equation (3.33) is used to give slower, more accurate convergence. 

3. Recursive Least Squares (RLS) Algorithm 

The recursive least squares, or RLS algorithm, is similar to the well 
known Kalman filter except that it is used to estimate model parameters rather 
than the system state. The tensor form of the normal equations is not suitable 
for RLS implementation. The elements of the correlation tensor must first be 
rearranged in a two dimensional matrix. This can be accomplished by replacing 
the tensor outer products appearing in equation (3.25c) with matrix Kronecker 
products. The Kronecker product of an n>m matrix, A= [a A A ], and a p>q 
matrix. B = |b M ], is an npxmq matrix given by [Ref. 30.31], 



a u B a 12 B . . . a lm B 

3.2jD ^22^ • • • 



A £B - 



(3.34) 



a n) B a r , ; B • • a Iim B 

where the symbol ;£ is used to denote the Kronecker product operation. 

In order to rewrite the normal equations. (3.25c). in this matrix form, the 
covariant tensor of system parameters, H A( . A . must be put into a vector form 
by a lexicographic reordering of the elements. The resulting parameter vector is 
denoted H. The normal equations become 



E x(k) c< 



x(k-N) 



x(k) £ • • • £,x(k 




H 



E-jy(k)x(k): x(k .\)j 

If an a-Miiiiption of ( rgodicit \ i- made then tin 
l>o replaced with time averages. 



(3.35) 



-talistical averages in (3.35) can 



lim — r 
K K 



K I 

£ 



k- 0 



X(k)X T (k) 



H 



Inn K £ l y(k)X(k) 

K-cx K k .„ 



(3.36) 
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where X(k) •- x(k) $<:••• 0 x(k N) . For computational purposes equation (3.36) 
is approximated as 

K— 1 K— i 

E X(k)X T (k) = £ y(k)X(k) (3.37) 

k = 0 k = 0 



As a further notational convenience, the following two definitions are made, 



X T (0) 

X T (1) 



and, 



*(K) 



y(0) 

y(i) 



y(K) 



(3.38) 



(3.39) 



Using these definitions the normal equations can be written in the 
compact form 



X h T X h H = X. t Y 



(3.40) 



This last from of the normal equations is precisely the same as the one used by 
Goodwin and Payne Ref. 32: p. 176 for the derivation of the RLS algorithm. 
The derivation involves the use of the matrix inversion lemma Ref. 33: p. 
247] which replaces a matrix inversion at each iterative step by a simple scalar 
division. It is this simplification which yields the efficiency of the RLS algorithm. 
The RLS equations aro (Ref. 32: pp. 176-177 . 

Hk-i H h - Qk-iI.v(K) X 1 (K - 1 )H h ) (3.41a) 



Qk i 



P k X(K- I) 

I * X r (K- 1 )P k X(K- 1) 



(3.41b) 
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Pk + i = iI-QK^.X T (K-l)P K 



(3.41c) 



= IX^Xk,,]- 1 (3.41d) 

Equation (3.41a) deserves some comment. The term in the parentheses 
is normally known as the innovation. It is a scalar which represents the new 
information gained in the latest measurement. The term X t (K+1)H k is an 
estimate of the current output given the new input measurement and the old 
estimate of the system parameters. The innovation is thus an error signal. The 
term in front of the bracket, Q K+1 , is a gain vector. It gives an indication of how 
much faith is being put in the new information. If the gain is small, then the 
new estimate will be essentially the old estimate. Conversely, if the gain is large 
then the new estimate will depend to a great degree on the new information. 
When the algorithm has converged, the gain will be close to zero. If the system 
parameters change for any reason, the algorithm will not detect the change since 
it is ignoring the new • information. In order to circumvent this problem a 
weighting can be applied to the data, so that new data is artificially emphasized. 
Exponential and rectangular windows have been successfully employed for this 
purpose. The time constant used in the case of the exponential window is often 
called a forgetting factor since it ensures that the algorithm’s memory is finite. 
Use of windows will not be discussed further in this dissertation. The interested 
reader should consult reference [Ref. 32: pp. 179-185]. 

4. Simulation Results 

In order to investigate the validity of the theoretical results, several 
computer simulations were performed. Two examples are presented, representing 
two different classes of systems. Many other sy^reim were also tested, however, 
tin result ^ presented are typical of Those obtained from all the simulations. The 
FORTRAN programs written allowed nonlinearities of up to fourth order, but 
were limited to systems involving only zero or one delay. 

The first example was chosen to correspond exactly to the model 
equation (3.17). The system was excited by white, zero mean, uniform noise. 
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The problem was to identify the system parameters from only input and output 
measurements. The unknown eovariant H A x tensor was chosen to be 



.2 


-.4 


.03 


-.7 


0.0 


.5 


.35 


.11 


.9 


0.0 


.01 


1.3 


-.33 


.7 


0.0 


.43 


.81 


-.05 


.4 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 



(3.42) 



Therefore, the output equation consists of 16 nonzero terms, ranging from 
2x ,0 )(k)x (0 >(k-l) up to .4x^(k)x^(k — 1). This model will be called System I in the 
sequel. 

The other type of nonlinearity tested was one that was known to require 
an infinite number of terms. The particular example chosen for this was the unit 
step function which simulates a saturation type of nonlinearity. It is a convenient 
choice since an analytical solution can be calculated. 

The unit step function has different H >o tensors depending on the 
order of the model chosen. This is a result of the chosen coordinates not being 
orthogonal (this will be clarified in Section D.) The second, third and fourth 
order (nonlinearity order) models are given by 
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o 


H x 


- 


- 0 


"'o 




2 


4 






.5 


.75 0 . 



(3.43a ) 







1 45 -35 


H ^ 






Vi 




2 32 32 




. 



1 .10625 0.0 - 1.09375 



11 X 




1 45 0 35 o| 


A n 




2 32 32 



|.5 1.40625 0.0 - 1.09375 0.0 



(3.43b! 



(3.43c) 
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Notice that the unit step function has no memory. The expansion contains only 
odd powers of x(k) and a constant term since the unit step function can be 
written as a constant plus an odd function (the signum function.) We will refer 
to these three models collectively as System II in the sequel. 

The direct solution of the normal equations was tested on these two 
models as was the LMS algorithm. To date the RLS algorithm has not been 
verified. 

a. Direct Solution of Normal Equations 

To verify the direct (matrix inversion) method of solution of the 
normal equations, a FORTRAN program was written which estimated the 

correlation tensors, and performed the required matrix inversion. The program 
was written so as to allow the number of points used to approximate the 

correlations to be varied. The maximum power of x(k) was also made to be 
adjustable so that the effect of over or under modelling could be studied. The 
final adjustable parameter was the magnitude of the uniform, zero mean, white 
noise that was used to excite the system. Adjusting this last parameter affects the 
range over which the resultant model is valid. In an actual application, 

something about the range of expected system inputs would have to be known in 
order to select a good value for this parameter. 

The results for System I are given in Table 3.1 for several 
combinations of the three variable parameters. The results are remarkably good 
even with as few as 30 input samples. Since no measurement noise was 

introduced this is perhaps not surprising. 

Overmodelling did not present any problems. The additional terms 
were identified by the algorithm to be essentially zero. This is evident in Table 
3.1c. Undermodelling did introduce some inaccuracy. The coefficients identified 
are significantly different from the ones in the unknown system. However, the 
coefficients identified should represent the best second order approximation to the 
system, which will not in general be the same as the second order coefficients 
contained in the third order model. 



58 



a. 30 data points were used 

The maximum power of x(k) was 3 
The noise was uniform on (-1,1) 
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b. 500 data points were used 

The maximum power of x(k) was 3 
The noise was uniform on (-10,10) 
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c. 500 data points were used 

The maximum power of x(k) was 4 
The noise was uniform on (-1,1) 
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d. 500 data points were used 

The maximum power of x(k) was 2 
The noise was uniform on (-1,1) 



Hi 



.20016 

.759151 

-.033444 



.80030 

1.5108 

1.6619 



.037553 

.069939 

.24457 



. TABLE 3.1: System I Resiilts Using Full Matrix Inversion. 
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The results for System II are shown in Table 3.2. Significantly more 
points were required to accurately identify System II than were needed for 
System I. The unit step function cannot be represented exactly by a finite 
number of polynomials so it is not surprising that the solution is not precise. 
There is another factor that contributes to the poor accuracy. The choice of 
functions used as coordinates, equation (3.14), leads to ill conditioned equations 
(see Strang [Ref. 34: p. 135].) This problem will be corrected in Section D. 
b. Simulation Using LMS Algorithm 

To verify the LMS algorithm a FORTRAN program was written, 
which used equation (3.30), to adaptively identify the covariant tensor H Aq Xfi - 
The program allowed the convergence parameter. //, and the magnitude of the 
uniform, white noise to be varied. The results of the simulations are presented in 
Table 3.3 and 3.4 for System I and II respectively. Equation (3.33) was used to 
bound the convergence parameter, . The input excitation noise was chosen to 
be uniform on (-1.1) resulting in a bound for the convergence parameter of 
0< /i <0.3558. 

In general convergence was slow. The linear and "close to linear" 
terms were identified most rapidly. The highly nonlinear terms, involving high 
powers of x(k) or x(k-l) and their cross terms, were last to be identified. Their 
accuracy never reached that of the lower order terms. The algorithm was very 
sensitive to the setting of the convergence parameter. A value smaller than 
the bound predicted above was used to achieve satisfactory performance. 

D. GENERALIZED COORDINATE SYSTEMS 

In Section B, a coordinate system was introduced that was closely related to 
the Yoherra series (equation (3.14).) ThL system wa> subsequently used in the 
remainder of Section B and in Section C. There* is realh little motivation for 
choosing this particular set of coordinates. In fact there are very compelling 
reasons to search for other sets of coordinates. The set (3.14) can lead to poorly 
conditioned sets of equations (see Strang Ref. 34: p. 135]). a fact that was 
mentioned in the last section. In higher order systems this can become a serious 
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a. 1000 data points were used. 

The maximum power of x(k) was 3. 
The input noise was uniform on (-1,1) 



(Hi 



.49278 1.4301 

-.037083 -.10437 

-.025155 -.072552 
.084649 .080036 



.000519 - 1.1389 

.067337 .18706 

-.035273 .10310 

-.13334 -.15251 



b. 15000 data points were used. 

The maximum power of x(k) was 3. 
The input noise was uniform on (-1,1). 



|H V J = 



.50380 1.4131 -.0093646 - 1.1029 

-.0064558 -.024451 .030395 .026762 

-.0016623 -.021414 .0044444 .030779 

.00015024 .024709 -.029871 -.018542 



c. 500 data points were used. 

The maximum power of x(k) was 2. 
The input noise was uniform on (-1,1). 



H, A 

■ <ri 



.47416 .77747 .00063628 
.0049366 -.0018654 0025596 
-020319 -.10698 .076991 



d. 5000 data points were used. 

The maximum power of x(k) was 2. 
The input noise was uniform on (-1.1). 



11 



AflA i 



.49480 .75558 .0025051 

.017591 .0011427 .024714 

(25590 0019598 052117 



TABLE 3.2: System II Results Esin" Full Matrix Inversion. 
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The maximum power of x(k) was 3. 

The noise was uniform on (-1,1). 

The convergence factor, n. was chosen to be .2. 



a. After 100 Iterations 



I H vJ = 



b. After 300 Iterations 



c. After 500 Iterations 



i H vJ = 

d. After 1000 Iterations 



e. After 1800 Iterations 



.23961 -.58084 

-.22910 .62263 

-.42404 .85335 

-.29008 .45359 



.21191 

.49735 

.11778 

.43573 



.19133 

.52472 

.032839 

.44209 



.18995 

.51438 

.028302 

.43507 



.20060 

.50360 

.017045 

.43879 



19456 

.48923 

0079652 

45295 



-.54901 

.55747 

1.3372 

.41638 



-.45457 

.65421 

1.2483 

.55546 



-.35146 

.55292 

1.2505 

.60299 



-.39481 

.46156 

1.2743 

64136 



-.38962 

.13032 

1.2786 

.67660 



-.089402 

.055704 

-.34090 

-.086633 



.088102 

.047074 

.28989 

082106 



.028912 
.10901 
-.32721 
- .070640 



.030956 

.11714 

-.33206 

.055581 



-.69479 

.71551 

.75055 

.53219 



-.68228 

66728 

.73386 

.56355 



70631 
.74664 
.75039 
.64 572 



.70562 

.75506 

.72930 

.61097 



TABLE 3.3: System I Results Using LMS Algorithm 



G2 



The Maximum power of x(k) was 3. 

The input noise was uniform on (-1,1). 

The convergence factor, was chosen to be .15. 



a. After 100 Iterations 



|H V J 



b. After 300 Iterations 



iH V ,l = 



c. After 500 Iterations 



!H 




d. After 1000 Iterations 



[H 



.27229 


.64437 


-.029974 


-.30266 


.021558 


-.13134 


-.16022 


-.10782 


.062328 


.049563 


-.069966 


-.17809 


.26624 - 


.00231392 -.093028 


-.048308 


.49886 


.96021 


-.11685 


-.88415 


-.010596 


.10195 


-.060408 


-.12604 


.036964 


.15373 


-.17940 


-.32627 


.14847 


.059296 


-.049165 


-.10986 


.54880 


1.3925 


.010187 


-.88497 


-.091216 


.14549 


.10978 


.0030069 


-.081292 


.23020 


-.14849 


-.33871 


.0091398 


.021484 


- 012741 


-071461 


.53253 


1.5226 


.044331 


1.0004 


.13698 


.060273 


1 5947 


.036425 


.16666 


.16941 


- .096854 


-.35439 


068181 


.018017 


.081221 


0062654 



c. After 1700 Iterations 



.48682 


1 2719 


- .058312 


.95098 


13368 


.025525 


.012328 


.037707 


.096072 


.20474 


-.031694 


-.21534 


05853 


.017710 


-.03364 5 


.025030 



TABLE 3.4: System II Results Using LMS Algorithm 
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problem. It is shown in this section that proper choices of coordinates lead to 
diagonal correlation matrices so that no costly matrix inversions are required to 
solve the normal equations. This makes the identification process almost trivial. 
It was stated in Section B.3 that the input observation vectors could be 
considered to describe a curve in a hypothetical (p+l)-dimensional space. The 
system output is estimated based on this curve (equation (3.17).) This chapter 
generalizes this idea to include cases where we have two curves, one depending on 
present and past inputs, the other depending on past outputs. Finally, the 
chapter concludes with several non-trivial numerical simulations. 

1. Choices of Coordinate Systems 

A point in (p+l)-dimensional space is defined by the variables 
x°, x 2 , ..., x p . In Section B.3 (equation (3.14)) these variables were chosen to be 
parametric functions of the variable x(k), the system input. The particular 
choice presented there was chosen to correspond to the Volterra series. It is 
desirable to pick a system of coordinates that ensure a diagonal correlation 
tejisor. as this allows solution of the normal equations without matrix inversion. 
A tensor T is diagonal if its components obey the rule 



pA 0 



/ 




< 



for A 0 = A . 



0 

V 



'b ~ ^ n+i y 

*> 

otherwise 




- 



(3.44) 



Note that only tensors of an even order (possessing an even number of indices) 
can possibly be diagonal. In general components satisfying the upper condition of 
equation (3.44) are called diagonal elements, or diagonal components. 
Components that are not diagonal are called off diagonal. 

Two conditions are required in order to ensure a diagonal correlation 
tensor. The first is a result of the following theorem, 
a. Theorem 3.1 

If a set of functions { f 0 (x). fj(x) f N (x)} are defined with the property 

that 
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K A for X — fi 
0 for X =* /i 



(3.45) 



f f>(x)f^(x)w(x)dx = 



where w(x) is a positive weighting function, then the set of random variables. Z\ 
defined as 



Z° = f 0 (X) 

Z 1 = f,(X) 

(3,46) 



V - f p (X) 

where X is a random variable with probability distribution 

p x (x) = Cw(x) (3-47) 

are uncorrelated. In equation (3.47) the constant. C. and weighting function. 
w(x), must be chosen so that p x (x) satisfies the definition of a probability density 
function. 



b. Proof 




e|z a z"| = e|c(x)c(X)| 


(3.48a) 


oc 

= Jb(x)f„(x)px(x)dx 
— 00 


(3.48b) 


oc 

f fjdx)f„(x)Cw(x)dx 

- OC 


(3.48c) 


00 

= C f (x) w (x)cix 

— OC 


(3.4Sdj 


Substiniting equation (3.45). yields 

■H ’ t« ; 7 


(3.49) 



Choosing a set of functions in accordance with equation (3.45) and 
ensuring that f 0 (x) is a constant function (ie: is a constant) is the first condition 
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that must be satisfied in order to obtain a diagonal correlation tensor. The 
system must be excited with samples drawn from a random process with 
probability distribution given by (3.47). In addition, if the random process is 
chosen to be zero mean and strictly white (ie, independent which implies 
uncorrelated), then the correlation tensor will be diagonal. This last condition is 



equivalent to requiring that 




E-^(m)x(n)j- = <5(m-n) 


(3.50a) 


p(x(n),x(m)) = p(x(n))p(x(m)) 


(3.50b) 


and 




o 

II 


(3.50c) 



where x(m) and x(n) are two input samples taken at times m and n respectively. 
It is a straight forward matter to show that if these conditions are met. the 
correlation tensor will be diagonal. 

The conditions presented above imply that different sets of 
coordinate functions should be used depending on the probability distribution of 
the noise used to excite the system. Different noise distributions have the effect 
of weighting the error differently. Consider that a Gaussian noise will contain 
samples of all amplitudes while uniform noise is bounded. If. for example, the 
system contains a saturation type of nonlinearity, the uniform noise may not 
detect its presence if its maximum amplitude i- hot sufficiently large. 
Theoretically. Gaussian noise contains samples of all amplitudes and will excite 
all modes. On the other hand it may be known that the -y-‘<*m input never 
exceed* a certain maximum value and ^o a bounded input will lx suitable. 

If two models of a system are constructed, in two different coordinate 
system* but using the same input noise (only one* set can possibly lead to a 
diagonal correlation tensor) then the two solutions will be equivalent. One 
solution can be transformed into the other by performing a change of coordinates 
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in accordance with equation (2.22). Therefore, it always makes sense to identify 
the system using the coordinate system which leads to a diagonal tensor. Other 
representations can then be calculated as desired. Solutions obtained by using 
different input noise density functions are not equivalent even if the coordinate 
systems used were the same. The effect of the different distributions is to weight 
the errors differently. A uniform input noise will weight the errors equally, while 
a Gaussian noise will emphasize the importance of errors made for small inputs. 
In general, transformations between solutions obtained using different excitation 
noise distributions cannot be found. Choosing an appropriate distribution 
requires knowledge of the expected system input signals. 

The Hermite polynomials lead to a diagonal correlation tensor if the 
input is white, zero mean, Gaussian noise. Similarly, Legendre polynomials 
should be used in the case of uniform noise. It is convenient to normalize the 
coordinate functions so that the diagonal components of the correlation tensor are 
all ones. To identify the system parameters, only the cross-correlations of the 
right-hand-side need to be calculated. This fact was first discovered by Lee and 
Schetzen [Ref. 10]. 

2. Recursive Models 

Recursive models have been used very successfully for modelling linear 
systems [Ref. 35]. Among their advantages is infinite memory, and the ability to 
model a system without knowledge of its input. The latter property allows these 
models to be employed in such areas as speech processing where input signals are 
difficult or impossible to measure. The assumption is made that the input is 
white noise. 

Recursive discrete-time nonlinear expansions have also been proposed 
(see for example !Ref. 14.15.16]). Recursive models poises infinite memory, and so 
may require fewer terms to accurately represent long memory systems. However, 
nonlinear recursive models also posses infinite nonlinearity. To understand why 
this is true, consider the following example. 
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a. Example 3.2 

Consider the following recursive, discrete, nonlinear system 

y(k) = ay 2 ( k — 1 ) u(k) (3.51) 

where u(k) = bS (k) (where £(k) is a unit sample) and where a and b are arbitrary 
constants. 

The output (y(k)) of the system for k = 1, K is 



y (o) = b 


(3.52a) 


y ( 1 ) = ab 2 


(3.52b) 


y ( 2) = a(ab 2 ) 2 = a s b 4 


(3.52c) 



y(K) = a 2 * -1 b 2 " (3.52d) 

It is clear that the nonlinearity of the system increases with time. Unlike a linear 
system, stability in a recursive nonlinear system is determined not only by the 
system parameter, a. but also by the input function. It is also difficult, in 
general, to predict for what classes of input a particular system will be stable. 
By analogy to the linear case we will refer to these types of models as Auto- 
Regressive or AR. 

It is possible to also expand the output as a combination of both past 
and present inputs and past outputs. This type of model wa* propo>ed by Parker 
and Perry f Ref. 13]. We will refer to this type of model as an Auto-Regressive. 
Moving- Average or ARMA model. It will have the *ame stability problem* a* 
does the AR model. 

Equation (3.17) can be used to model an AR nonlinear system* if the 
proper choice is made for the observation vectors. Using t lie >ame coordinate 
functions as given in equation (3.14) but using y(k - i), i = 1,....N, as an input 
parameter, yields appropriate observation vectors. That is. 
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y A '(k-i) - y(k-i)j tAi) (3.53) 

defines the components of the i-th observation vector. The model equation, 
equivalent to equation (3.17), becomes 

y(k) = y A ‘(k-l)...y AN (k-N)H Ai . . , N (3.54) 

where H A ... A is again a (p+l)-th order covariant tensor containing the system 
parameters. This model is an extension of the familiar linear autoregressive, or 
AR model. 

The normal equations derived for the nonrecursive case can be used 
to solve for the H A] . An tensor in this case as well. The identification process is 

described in Figure 3.5. The assumption is made that the system is recursive and 
driven by a white noise, u(k). Its output can then be described by 

y(k)= y A, (k-l)..y AN (k-N)H Ai . AN + u(k) (3.55) 

This output signal is delayed and fed into the analysis model. The error signal 
e(k) is given by 

e(k) = y(k) - y(k) (3.56a) 

= y Al (k-l)...y AN (k-N)H v - u(k) - y A ‘(k ])...y A "(k-N)H Ai . . An 



(3.56b) 

When the model parameters exactly equal the actual system 
parameters the error signal will equal the input white noise. For this rea>on the 
analysis model is often called a whitening or bleaching filter. The analysis 
model is nonrecursive. It uses past values of the system output to make a 
prediction of the present system output. The normal equations (3.25c) apply to 
this situation with the observation vectors, x(k-i). replaced by the vectors defined 
in equation (3.53). The normal equations for the recursive nonlinear model can 
be written as 
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6 




£ 



Figure 3.5: Recursive Model Identification 
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(3.57) 



Ely l (k-l)...y N (k-N)y'* 1 (k- " N (k-N)fH Ai 



= E{y(k)y Ml (k~l)...y" N (k-N) 



These equations are a tensor equivalent of the Yule- Walker equations. The 
corresponding mean-square error is 




- Ely (k)y *(k— 1) 



y" N (k- 



N) H; 



(3.58) 



In this case it is not obvious that one set of coordinates will yield 
better results than another. The correlation tensor appearing on the left-hand- 
side of equation (3.57) cannot be guaranteed to be diagonal since the probability 
distribution of y(k) cannot be controlled. Techniques must be devised which 
choose the coordinate system "on line" as the statistics of y(k) are determined. 
For example, in the linear lattice the coordinate system is chosen by 
orthogonalizing the sequence y(k) using a Gram-Schmidt procedure. In this way 
the coordinate system is not determined until run time. Extension of these ideas 
is left until Chapter 4. 

The tensor model presented is also' suitable to represent a nonlinear 
ARMA model. This type of model takes into account all available information 
(input and output) and so should be more accurate. It may also lead, in some 
cases, to a lower order solution than either the AR or MA model. An ARMA 
tensor model can be written as 

y(k) = x A# (k)...x A “(k-M)y'‘ , (k-l)...y ,,N (k-N)H A# . ^ „ N (3.59) 

Relations analogous to (3.57) and (3.58) for the normal equations 
and the minimum mean-square error for the ARMA model can be obtained in a 
straight forward manner. 

3. Simulation Results 

The concepts developed in Section D of this chapter were verified using 
FORTRAN simulations. The coordinate functions were chosen to be the 
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normalized Legendre polynomials because of the ease in generating uniform noise 
on a digital computer. By assuming that the correlation matrix was diagonal, 
the solution to the normal equations was determined without matrix inversion. 
This solution was verified by generating the correlation on the left-hand-side of 
the normal equations and performing the required matrix inversion. The LMS 
adaptive algorithm was also tested using the Legendre polynomials. 

The Legendre polynomials used are given by 



f 0 (x) = 1.0 


(3.60a) 


X 

(eo 

/> 

11 


(3.60b) 


L(x) = (x« - |) 

V X u 


(3.60c) 




(3.60d) 



These functions are normalized so that the correlation tensor will have ones along 
the diagonal when excited with zero mean white noise which is uniform on (-1,1). 

The simulation models used were similar to those used in Section C. The 
coefficients for System I were identical to those given in (3.42). although thi- 
time they are coefficients of Legendre polynomials so they do not represent the 
same system. System II was again a unit step function which has a tensor 
representation, in terms of the Legendre polynomials given by [Ref. 26: p. 1526], 



H a 




1 n'S Q _ x 175 1 






2 4 80 



(3.61a) 



o.. r > 0 . 4 SX 0 1 3 o.o o.i6r»:;r>9 



(3.61b) 



The results of the simulations for System I are presented in Tables 3.5 
and 3.6. Results for System II are given in Tables 3.7 and 3.8. In all cases it i> 
obvious that the simulation results are very good. 
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15000 sample data points were used. 

The maximum power of x(k) was 3. 

a. Solution Without Using Matrix Inversion 

.19158 -.38475 .0041187 

.49067 .38457 .089293 

-.012253 1.2925 -.35452 

.41049 .85202 -.074247 




b. Solution Using Matrix Inversion 



! H v.l 



.20000 -.40000 .030000 

.50000 .35000 .11000 

.01000 1.3000 -.33000 

.4000 .81000 -.050000 



-.67056 

.93948 

.71923 

.43011 



-.70000 

.90000 

.70000 

.40000 



TABLE 3.5: System I Model Parameters Obtained 
a. Using No Matrix Inversion, and 
b. Using Matrix Inversion. 



The maximum power of x(k) was 3. 

The convergence factor was chosen to be .01. 



a. After 200 Iterations 





.14678 


-.45852 


-.036398 


-.72277 




.42720 


.25876 


.018320 


.83079 


(H V J = 


-.077103 


1.1898 


-.45388 


.61647 




.34264 


.68052 


— .16955 


.30546 


b. After 500 Iterations 




.20013 


-.40047 


.030509 


-.70066 


II 

X 


.50038 


.35077 


.11136 


.90025 


.0099732 


1.2996 


-.32958 


.70026 




.43008 


.81014 


- .050063 


.39875 


c. After 1000 Iterations 




.20000 


.40000 


.030007 


-.70000 




.50001 


35000 


.11001 


.8999 


a = 

''O'M 


.010008 


1 3000 


32999 


.69999 




.43001 


.81000 


- .049986 


.40000 



TABLE 3.6: System I Model Parameters Using LMS Algorithma. 
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a. 500 data points were used. 

The maximum power of x(k) was 3. 



No Matrix Inversion Solution 



[H Vl ] = 



.52104 

-.011098 

.0011160 

.040052 



.43930 

-.0046220 

-.010606 

-.0071009 



.0090070 

.014464 

-.010466 

.058263 



-.17515 

.020122 

.021310 

.052566 



Using Matrix Inversion 



:»v,i = 



.50157 

-.0084508 

.0017041 

-.0066652 



.43054 

.0032387 

-.0009167 

-.0024676 



-.011153 
- .0000970 
.0062532 
.012192 



-.16530 

-.0006363 

-.0059126 

.0014058 



b. 15000 data point were used. 

The maximum power of x(k) was 3. 



No Matrix Inversion Solution 




.50050 

.0015150 

-.0071103 

.0022021 



.43363 

.0003272 

-.0004279 

.0018209 



-.0043279 

-.0012861 

.012933 

-.0021332 



-.16772 

-.0003738 

.014358 

-.0060674 



Using Matrix Inversion 



H 



A .A . 



.50089 

0018213 

0036101 

.0003990 



.43306 

-.0005741 

.0008143 

.0001746 



-.0001744 

- .0014410 
. 004381 1 

- .0003974 



.16722 

.0000789 

.0008015 

.0006467 



TABLE 3.7: System II Model Parameters Obtained 
a. Using No Matrix Inversion, and 
1). Using Matrix Inversion. 
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The maximum power of x(k) was 3. 

The convergence factor, n, was chosen to be .005. 



a. After 300 Iterations 



( H A 0 aJ - 



.46464 
.011790 
— .010943 
.0014149 



.41356 

.012146 

.014180 

-.0065302 



.0066999 

.0067764 

.020212 

.0093879 



-.16385 

-.0093346 

.018152 

-.0047952 



b. After 500 Iterations 



1«a. 



.50556 

-.015344 

-.0082813 

.012877 



.43020 -.00098824 -.16796 

-.010684 .020836 .0098449 

-.0052115 .0090727 -.010251 

-.0094530 -.0046525 .0072316 



c. After 1000 Iterations 



H, 



1 = 



.49868 .44629 .0053940 

.020577 .0091397 -.010215 

012857 .020072 -.0064996 

-.0020197 .0031219 -.0064275 



-.17370 

.0065185 

-.010541 

-.0037119 



c. After 1800 Iterations 



Hi 



.50115 

0000129 

0005666 

-0075634 



.43727 

-.0063172 

.0084999 

.0058037 



-.0059911 

-.012713 

-.0068103 

.15156 



-.16508 
.010275 
- .010135 
. 01*235 



TABLE 3.8: System II Model Parameters T'sing LMS Algorithm 
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IV. REVIEW OF LATTICE FILTER STRUCTURES 



This chapter reviews lattice Filter theory. These Filters were First proposed in 
connection with the linear prediction of speech by Itakura and Saito [Ref. 36] in 
1971. They developed a new approach based on a partial correlation (PARCOR) 
coefFicient. Since that time the Filter structure they proposed has come to be 
known as a Lattice (or sometimes also called Ladder) Filter. The properties of the 
Filter have been exhaustively studied by many researchers [Ref. 37], Lattice 
filters have been successfully applied to problems in various disciplines. 

The great interest in the lattice approach stems from it's property of 
orthogonality. This property allows the Filter to be updated in order, without 
recalculation of all the previous, lower order, Filter coefFicients. Orthogonality also 
leads to a nice physical structure, a cascade of First order sections, and so is 
appropriate for efficient hardware or software implementations. Finally, it has 
been shown that the lattice owes its robust numerical behaviour to this property 
of orthogonality [Ref. 38: pp. 128-1361. 

This chapter begins with a derivation of the one-dimensional (1-D) lattice. 
The approach taken in the derivation is somewhat novel in that it begins by 
expressing the linear prediction in terms of an uncorrelated error sequence. This 
is regarded as a change of coordinate systems and used to develop the Levinson 
algorithm and the lattice filter. In the following section generalized forms of the 
Levinson algorithm and lattice are derived. It is shown that these also 
correspond to coordinate transformation. Next, the Schur algorithm, which is a 
method for generating the required filter coefFicients (Lattice parameters) 
directly, given only a knowledge of the correlation matrix, is reviewed. In the 
following chapter the generalized lattice filter is used to develop new. 
multidimensional (specifically 2-D) lattice filters. In Chapter 6. a new nonlinear 
lattice. filter, based upon this generalized lattice formulation, is presented. 
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A. 1-D LINEAR AUTOREGRESSIVE LATTICE FILTER 



We define the N-th order forward error sequence of an autoregressive model 
as 



e N (k) = y(k) - £ h A N y(k - A) (4.1) 

A= 1 

This can equivalently be written in tensor notation as, 

e N (k) - hjV (4.2) 

where 

[y A j T = |y(k) y(k-l) ••• y(k-K)) (4.3) 

and 



[h A N ] = [1 -h, N -h 2 N • • • -h N N 0...0] (4.4a) 

where for convenience, we have made all vectors of length K + l (ie; A = 0. .... K), 
where K is some arbitrary maximum length. Note that, 

h 0 N = 1 * (4.4b) 

The y ; can be considered to comprise a coordinate system in an K + l 
dimensional space. The vector y A ; represents a single realization from a random 
process. As mentioned in the Chapter I, there will in general be many such 
vectors corresponding to all the possible realizations of the random process. 

Because the error. e N (k). is a scalar (invariant) it must remain unaltered 
regardless of the choice of coordinate system. We may. therefore, write 

, s (k) = K a V (4.5) 

where the y A are defined as 
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y( k ) 

y(k-l) 

y(k-2) - hoVfk-l) 
y ( k 3 ) - h 2 y(k — 2)h 0 2 y(k- 1) 



y(k-K) -h^Mk-K+l)- ••• - h 0 K_1 y(k-l) 



and 



(4.6) 



[K A N 'j = (1 -K, -K 2 • • -K n 0...0] 



(4.7) 



where the h’s are chosen so that the components, y A , are uncorrelated. The 
stochastic form of the Gram-Schmidt procedure can be used. A discussion of this 
method can be found in [Ref. 39: pp. 382-383]. By uncorrelated we mean 



eVV 



fb x ' for A ' - /x 
\0 otherwise 



(4.8) 



The reader familiar with lattice structures will recognize the components of y A as 
the backward prediction errors. That is they are the errors in predicting y(k-X) 
from the next N-l values: y(k-N-f-l). ..., v(k-l). 

It is a straight forward matter to solve the prediction problem given by 
equation (4.5) (ie; solve for the K’s) because of the uncorrelatedness of the chosen 
coordinate system. Using an approach similar to that presented in Chapter III. 
Section C. a set of normal equations can be formulated for this problem. In this 
case, however, the correlation matrix is diagonal so that there is no inversion 
necessary to obtain the solution. Minimizing the mean square value of the error. 
e N (k). given by (4.5). with respect to the K A N s. yields 



K 



\ 



i:jy(k)\ l j 

u {(> 1 rj 




(4.9) 



Having obtained a solution to the orthogonal problem, we wonder if it cannot 
be employed to advantage to simplify the calculations required to solve the 
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original autoregressive problem. From equation (2.13). we know that the 
relationship between the old and new components can be expressed as 



dy 



A ' 



dy A 



where 



8y A 

9y A 



1 0 0 

0 1 0 

o -ho 1 1 

0 -h 0 2 -h, 2 



0 0 
0 0 
0 0 
0 0 



o -ho K " 2 -h,*"* • • • 1 0 

0 -h 0 K 1 -h^- 1 • • • -h^ 1 1 



Now. since 



e N (k) = yV = y A 'K A N - 



(4.10) 



(4.11) 



\ 



(4.12a) 



= V 



ay ; 



K a N ' 



(4.12b) 



then 






dy A 

0 V A 



(4.13) 



Equation (4.13) gives the relationship between the autoregressive model 
parameters and the K A N ’s. This result will be used in the proof of Levinson’s 
algorithm. 

1. Levinson-Durbin Algorithm 

In 1947 in his classic paper [Ref. 40 Levinson develops! a method for 
recursively solving the normal equations. Beginning with a zero order solution 
successively higher order solutions are calculated. This algorithm can be used to 
exploit the Toeplitz nature of correlation matrices of stationary random processes 
in order to reduce the required number of computations. The algorithm as 
presented in this work is actually a simplified version of Levinson’s original from. 
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It assumes that the equations being solved are the Yule-Walker equations ('Ref. 
41], see also Chapter 3. Section D) so that the right-hand-side of the equations 
contains only terms that will be present in the correlation matrix (left-hand-side) 
of the next higher order model. This simplification is due to Durbin [Ref. 42]. 
We will refer to this algorithm often simply as the Levinson algorithm, however, 
Durbin’s contribution is acknowledged. 

a. Theorem 4.1 Levinson’s Algorithm (Durbin’s Form) 

The autoregressive model parameters may be calculated recursively 
from the relation 

!h A N+1 ] = [h A N ] - K N .,S'h A N ] (4.14a) 

where 

!Ha N ] - [-h 0 N -h, N • • • -h,?., 1 0 • • • 0] (4.14b) 

and 

S[h A N ] = [0 -h 0 N -h, N • • • -h N N _, 1 0 • • • 0] (4.14c) 

The operator S has the effect of shifting the components. h A s one position to the 
right. Note that 

hj? - l 

For stationary processes the following simplification applies 

h A N - h N N _ A _, for A - 0 N (4.14d) 

b. Proof (by Induction) 

Using equation (4.13) for the first order model we have 



V h, 1 0 • ■ • 


0 


(4.15a) 


1 K, 0 


0 


(4.15b) 


] 0 ■ • • 0 


K, 0 1 0 • • • 0 


(4.15c) 


h A ° - KjS|h A °. 




(4 . 15d) 
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so that equation (4.14a) holds for N = l. We assume the recursion holds for the 
N-th order model 



!hf = !h A N " 1 ] - K N S|h A N-1 
From equation (4.13) we have 

H> n = K a n - A, A' = 0, K 

dy 

Therefore, 



(4.16) 



(4.17) 



-K, + K/ho 1 + K s h 0 2 +•••+ K N h 0 N-1 
-K 2 + Kjh, 2 + K 4 h, s +...+ K N hj N_2 



h A N = 



K n _ 

K n 



kAI 



N n N-2 



(4.18) 



0 



0 



Now. for the (N + l)-st order model 



1 

K, - Koho 1 * K s h 0 2 K N ^h 0 N 

Ko K 3 h { - K 4 hi^ K\4_|hj^ 



lT 



h 



X 



K\ — , 

~ Kn - 1 
0 



(4.19a) 



0 
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1 

-Kj K/ho 1 + Kjho 2 + ...+ K N h 0 N -' 


T 


^ o ^ 

O ^ 

1 


-K 2 + K,h, 2 + K 4 h, s +...+ K N h, N_1 




-hi 


-k n 

0 


“ K N + 1 


T N 

_h N-l 

1 

0 


0 




0 



(4.19b) 



\ 



- h A N - K N+1 S[h A N ] 



(4.19c) 



And so the desired result has been confirmed. 

The proof presented does not require that the process be stationary 
and so the condition of (4.14c) is not necessary. One is then faced with the 
problem of how to determine the h A N, s. This question is answered in the next 
section where a generalized form of the algorithm is derived. We note that the 
condition of (4.14c) implies the following 

h\ = K n = h 0 K (4.20) 

or that the second column of the matrix given in equation (4.11) contains all the 
lattice coefficients. 
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The following significant observation about equation (4.14a) is made. 
The transpose of the term within the second bracket on the right hand side 
corresponds exactly with the rows of the coordinate transformation matrix of 
equation (4.11). The Levinson algorithm, therefore, represents a recursive 
method for finding an orthogonalizing coordinate transformation. 

2. The 1-D Lattice Structure 

It has been shown, in the previous section, that the autoregressive model 
parameters can be calculated in a recursive manner using the Levinson algorithm. 
It was also shown that a model could be built in an orthogonal coordinate system 
where the model parameters were given by the K^’s. The question arises, can a 
filter structure be devised which represents y(k) in the orthogonal coordinate 
system? The answer is af^lrmative^ It will be shown in this section that the 
Lattice filter is the required structure for the stationary case. A more general 
solution is presented in Section B of this chapter. 

The desired result is obtained in a straight forward manner by 
multiplying both sides of equation (4.14a) by y\ This yields 

e N+1 (k) - e N (k) - K N+1 r N (k — N — 1) (4.21) 

where the quantity r N (k-N-l) = y A evaluated at A'=N+1. As was mentioned 
earlier in this chapter, the quantity r N (k-N-l), is generally known as the 
backward prediction error since it corresponds to the prediction of the point v(k- 

N-l) from the X future points y(k-l) y(k-N). 

Assuming stationaritv. we can use (4.14c) and (4.14a) to obtain 

-- Slh*' K n ., hr (4.22) 

which leads directly to the equation 

r N>1 (k — N— 1) - r N (k-N 1) - K N + 1 e* N (k) (4.23) 

Equations (4.21) and (4.23) define the Lattice form of the whitening filter (see 
Chapter 3, Section D). This is also sometimes referred to as the analysis model. 
The structure is illustrated in Figure 4.1. 
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In order to develop a synthesis model we need only to solve equation 
(4.21) for e N (k) since we know that e°(k) is equivalent to v(k). Therefore, the 
synthesis equations are 

e N (k) = e N+1 (k) + K N+1 r N (k— N— 1) (4.24a) 

r N+1 (k— N— 1) = r N (k— N— 1) - K N+1 e N (k) (4.24b) 

The resulting structure is as illustrated in Figure 4.2. Note that in this model, 
y(k) is indeed being expressed as a weighted sum of the backwards errors, which 
represent an orthogonal coordinate system. 

This concludes the discussion of the classical Lattice formulation. 

B. GENERALIZED ORDER UPDATE RECURSIONS 

In this section a more general linear prediction problem is considered. No 
assumptions are made as to the origin of the data. In fact, the data need not 
represent a time series at all, and certainly shift invariance is not required. The 
ordering of the data is simply chosen in some convenient fashion. The generality 
of this formulation will allow its application to multidimensional and nonlinear 
problems. 

1. Definitions and Formulation 

In this section quantities are defined that will be required to complete 
the statements and proofs contained in the remainder of the chapter. 

A realization of the random process Y is given by the vector 
y*;, 0 $ A $ K. 

The error. e k N l,. in predicting the element y k ~' from the previous N 
elements of Y is given by 

ekli - h; N (k-l)y A for A - 0 K (4.25) 

where 

h A N '(k^l) = iO • • • 0 hJi N _, -h k N : N%2 • • • -h k N 1 0 ■ • O' (4.26) 



85 



e° <k) 



e 1 (k) 



e 2 (k) 




Figure 4.1: 1-D Lattice Analysis Model. 



e° (k) e 1 <k) e 2 (k) 




Figure 4.2: 1-D Lattice Synthesis Model. 
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The norm of this prediction error is given by 



P N 1 

e k+l I 



11/2 



Ki, 



N (2 



(4,27) 



A normalized version of the forward error is given by 



- 

^k+ 1 



e k-f 1 

"n ' | 

e k + l I 



(4.28a) 



= a A N (k+l)y^ 



(4.28b) 



where 

a A N (k+l) = h A N (k-l) (4.29) 

e k + l 

The backwards prediction error, rjl N , is the error associated with the 
prediction of y k_N given the next N elements of Y. It is given by, 

r k N N = H A N (k-N)y A (4.30) 

where 

[h A N (k-N)] = [0 • ■ • 0 1 — hk-N^j -h k N i N+2 ‘ ’ ' ~h k N 0 • ■ • 0 ("1-31) 

Again, a normalized version can be defined as 

bt-N = ^ (4.32a) 



- b A N (k-N)y A 



where 



bi N - 



h A (k-N) 



-N 
r k — N 



and 



! I r k- N 




(4.32b) 



(4.33) 



(4.34) 
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In order to generalize the results of Section A of this chapter we need to 
introduce two different families of coordinate systems. We will refer to members 
of these families as either forward or backward local coordinates. The term local 
is used because a different coordinate system will be associated with every value 
of k and N. We define the (k+1) - indexed, N-th order forward coordinate 
system as 



V k-N+I lN- 1 k-N-t-2 

y n k-N-2> 



t N— l k 

hk y 



y A (k+l,N)j = 



(4.35) 




hkV 



The corresponding coordinate transformation from the unprimed coordinate 
system to this local forward coordinate system is given by 



dy A (k~ l.N) 
8 y A 
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1 1 




0 






l_o 


1 1 


I N — 1 

“ n k-N-2 


r N— i 

n k-N + 3 


L N-l 
- h k _ j 


I N -1 
^k 


1 


1 0 


1 


IN-2 
~ n k-N + 3 


IN-2 
“ h k- 1 


I N-2 
“ ^k 


1 


1 0 


0 


1 


T N — 3 
“^k-l 


v-* 


1 

i 


0 1 ' 




* 






i 

■ 0 


1 0 


0 


0 


1 


-W 


1 

i 


1 0 


0 


0 


0 


1 


1 


0 1 




0 






1 1 



( 4 . 36 ) 



where the O’s are zero matrices and the I’s represent identity matrices. 

Similarly, we define the (k-N)-indexed, N-th order backward coordinate 
system as 



,.o 



> -r. 



,k-N -1 



(k- 



,k-N+2 



_ ^k-N+iy 



k-N + 1 



( 4 . 37 ) 



- hk-i y 



N-r.k-i 



u N- 1 „ 



k-N + 1 



The corresponding coordinate transformation is given by 

by"' (k X.N) 
by* 
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it 






0 






l_o 


1 


1 


0 


0 


0 


0 


1 


1 


“W-N + l 


1 


0 


0 


0 


1 


1 


V, 2 

“ n k-N + l 


“hk-N+2 


1 


0 


0 


1 


0, 

1 

1 


v N-2 
~ n k-N+l 


i N-2 
“ n k-N+2 


i N "2 

~ n k — N + S 


1 


0 


0 

I 

1 

1 


1 


i N-l 
^k — N + l 


i N-l 
“ n k-N+2 


l N-l 
“ h k-N + 3 


V, N-l 
~hk-i 




1 


0 1 






0 






1 1 



( 4 . 38 ) 



These definitions are sufficient to state and prove the theorems presented 
in the remainder of this chapter. 

2. Generalized Levinson Algorithm 

The Levinson algorithm (equation (4.14)) can be extended to recursively 
compute the forward and backward prediction coefficients defined in equations 
(4.26) and (4.31). In this section two forms of the algorithm are presented and 
proved. First a non-normalized version is introduced, then using this result a 
normalized algorithm is developed. 

a. Theorem 4.2: Generalized Levinson Recursion (Regular Form) 

The forward and backward prediction coefficients defined in 
equations (4.26) and (4.31) can be updated using the following recursive 
equations. 



hK->(k-l); = h A N (k+l)' - p^hjftk-N); 



r k-N 



(4.39a) 



'hi '* 1 (k- N)j = hi N (k-.\) - 



'n-> h A N (k-l) 



**k-\ 



e k N rl 



where 



P\ + i - 




(4.39b) 



(4.40) 



b. Proof 

The forward and backward prediction errors are scalars so their 
representation is identical in all frames of reference. Therefore, we can write 
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e k+i = h A N (k+l)y A 



(4.41a) 



= K^(k-l)y A '(k+l,N) 



(4.41b) 



and 



r„ N _ N - h A N (k-N)y A 



= K/b(k-N)y A "(k-N,N) 



where 



K A N -(k+l) = h A N (k+l)- 



9v a 



9y A '(k+l,N) 



= 0 • • • 0 -K k N _ N ^ • • • -K k N 1 0 • • • 0] 



hfik+i) = K ) -:(k4-i) ay> >t 1 J i ) 

9y A 



and 



K A N -(k-N) = h A N (k-N)— - 



9y A 



9y A (k-N.N) 



= |0 • 0 1 -K k 1 N _ 



- k k N 0 • • • O' 



h A K (k-N) = K A N -(k-N) 8yA lk A N ’ N ) 

0 v 



(4.42a) 

(4.42b) 

(4.43a) 

(4.43b) 

(4.44) 



(4.45a) 

(4.45b) 

(4.46) 



The normal equations in the primed and unprimed coordinate 
systems can be solved for K A N -(k-l) and K A N b(k-X). This is once again (see Section 
A of this chapter) a straight forward matter because the correlation matrices. 
[E{y A (k-l.N)y M (k~l,N)}] and [E{y A (k X.XJy'' (k N.N)}]. will contain only 
diagonal elements. The solutions are 



K A (k-f 1] 
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Ejy(k- l)y k-N ~ 1 (k-t-l,N)| E'Sy(k-t-l)y k (k — 1,N)| 

0 • • • 0 7 r— - 7 7 1 0 • ' • 0 



E<(y k - N+ 1 (k+l.N )) 2 



Ey(y k (k— l.N )) 5 



KA N "(k-N) = 



(4.47) 



E*|y (k — N)y k_N+ 1 (k- N,N) f E^y (k- N)y k (k-N,N) > 

0 • • • 0 1 7 r— ^ - —7 r — — 0 • ' • 0 



E^(y k - N ^ 1 (k-N,N )) 2 



E)(y k (k — N,N )) 2 



(4.48) 



Consider the (k-N)-th term of the (N+l) order version of (4.47). 



Kk-V(k-l) = 



Eiy(k+ l)y k- ' (k+l,N) 



l|(y k - N (k-l,N)) 2 j 



(4.49a) 



E'i v ( k - 1 ) r kN ! 



E<!(r k N : N ) 2 



(4.49b) 



r- I N N 1 

kSe k ~ir k n ‘ 



/ 



E*i ( r k - % ) 



(4.49c) 



N 

e k~l k 

— P N*l 

r k N 



(4.49d) 
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Similarly, 



Kk N + V(k-N) = ' '- v -— 

e k+l I 



(4.50) 



For the N = 1 model the forward error prediction coefficients are 
(from equation (4.44)) 



[h A x (k+l)] = Ki-(k+l) By ' 

dy 



(4.51a) 



= [0 • • • 0 -K k ! 1 0 • • • 0] 



(4.51b) 



= io • • • 0 



Equation (4.39a) yields 



ek ° +l! 1 - Pl k 1 0 • • • 0] 



I ! r k° I I 



[hi(k+l)] - [h A °(k+l)] - /> 1 k |h A °(k)| 



e k°+i 



r k 



(4.51c) 



(4.52a) 



= [0 • • • 0 1 0 • • • 0 ] - pi 



0 0 
e k-rl 



r k 



[o • • • 0 1 o • • • 0 



(4.52b) 



= i° ••• 0 i o 



(4.52c) 



Therefore, the recursion of equation (4.39a) holds for N=l. We assume it is valid 
for N and verify it for N+l. From equation (4.44) we have 

9y r (k- l.N) 



h A N -‘(k-l). = Kr‘(k,l) 



N - 1 / 



d V' 



(4.53a) 
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0 



T 



0 

-Kk-N^k+l) 

— KjJL'5i 1 + i(k+l) + hk-N+iKk-N^k+l) 



-K k N+1 (k+l) + hX^k+l) + ••• + hW.tftk+l) 

1 

0 



(4,53b) 



= [hi"(k^l) - K k N :V(k-l) 0 • 0 1 -h k - N+1 • • • -h k N 0 • ■ • 0] (4.53c) 

N 

= h A N (k- 1) e V* /»i„'h A N (k-N) (4.53d) 

F k \ 

This completes the proof of equation (4.39a). Using similar 
arguments the backwards recursion of equation (4.39b) can be verified. 

c. Theorem 4.3: Generalized Levinson Recursion (Normalized Form) 

The normalized prediction error coefficients defined in equations 
(4.30) and (4.34) can be recursively updated using the following recurence 
relations: 



a^Mk-l) 
b N ^(k N) 






«A S (1<-1) 

fcJ(k-N) 



w here 



e(PN k -.) : 

V 1 ( P N - 1 ) ‘ 



1 k 

1 - l 

k 1 

P\ - l 1 



and px+i is given by equation (4.40). 



(4.54) 



(4.55) 
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d. Proof 



The proof follows directly from the non- normalized version of the 
generalized Levinson algorithm (equations (4.39)) and the definitions of the 
normalized prediction error coefficients (equations (4.29) and (4.33)). 

Using (4.29) in (4.39a) we obtain 

a A N+1 (k+l) 





a A N (k+l) 




r£ N 1 | b A N (k-N)] 



(4.56a) 



e k N +l 



i „ N + l 
! ; e k-i 



- a A N (k+l) - PN +1 b A N (k-N)] 



Similarly, using (4.33) in (4.39b) we obtain 



br’(k-N) - 



r k-N 



_ N + 1 
r k-N 



■|b A ft (k-N) - , N * +1 a A "(k-rl): 



e k + l ! 



The proof of the fact that 

I I r k % | | i 



e 



N + l 
k-1 



! j N + l 
r k- N 



\A ~ (p n-i ' 2 



(4.56b) 



(4.57) 



(4.58) 



is relegated to the next section. 

The initial values for forward and backward normalized 
autoregressive parameters are obtained by setting X=0 in equations (4.29) and 
(4.33). This yields 



a A °(k+l) =;()•■ 


• 0 - 

i 


i 

; y k "‘ 


— 0 • 


• ■ o; 


(4.59a) 


b A °(k) = 0 • ■ 


■ • o - 


i 

i _.k 


- o • • 


■ 0 


(4.59b) 



We note the similarities in the two generalized forms of the Levinson 
algorithm presented in this section with that presented earlier in this chapter. A 
little reflection will convince that equations (4.14) are simply a special case of 
equations (4.39). 
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3. Error Order Update Recursions 
a. Theorem 4.4 



The (N+l) order errors can be calculated from the N-th order errors 
through the recursion 



-N + l 




— N 


e k + l 


- ®(Pn + i) 


e k + l 


it 


f n 

r k-N 



where 0 (/>n+i) is defined in equation (4.55). 
b. Proof (Outline) 

Using the normalized form of the Levinson algorithm (equation 
(4.54)) and the following relationships 

ejti = a A N (k+ 1)/ (4.28b) 

e k N .V = a A N+1 (k+l)y A 



F k N - N = b A N (k-N)y A 

~ b A N + 1 (k-N)y A 

equation (4.60) is easily verified. 

We now return to the proof of equation (4.58). That is. 

e k N ,i' 1 _ r£ K _ 1 

! " r£y, " s/iq^Tjl 

Working with the term on the left-hand side of the equation 



(4.32b) 



(4.58) 





N - 1 \ 2 



(4.61a) 



k - 1 N ) 



E.v^V-! 



Ey k+1 e k 1V 



T7T 



(4.61b) 
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1 



THT 



„ k pJ k-^l-N 
+ r k— N 



e k-l I 



1 - 




(4.61c) 



1 

\/l _ (pn+i ) 2 



(4.61d) 



Similar arguments can be used to verify the other relationship given in equation 
(J.58). 

Figure 4.3 symbolizes a single stage of the recursions of equation 
(4.60) while Figure 4.4 illustrates a third order analysis model. 

The interesting feature of equations (4.60) is that they do not make 
any assumptions about the nature of the given data. The data values need not 
be delayed versions of each other, as is the case for the autoregressive model. 
Any set of data values can be used. This fact will be significant when we deal 
with 2-D and nonlinear lattices. 

4. The Generalized Schur Algorithm 

In this section an algorithm will be presented which allows the 
calculation of the partial correlation coefficients in a direct manner. It will be 
shown that knowledge of the correlation matrix is sufficient to calculate all the 
reflection factors and thus solve the normal equations (by use of tin Levinson 
algorithm). The method used has come to be known as the Schur algorithm Ref. 

31 - 

In order to obtain the desired result we must introduce two new 
quantities, defined as 
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r 



N+l 

k-N 



\ 
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Figure 4.3: Generalized Lattice Filter Sections 
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Figure 4.4: 3-rd Order Generalized Lattice Filter 
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Q N-t-i(k+l) = Ek N ; iy M = a*(k-l)R" A 



( 4 . 62 ) 



^N+i(k— N) = Ek N _ N y4 = b, N (k-N)R^ 



where 



(4.63) 



R" A - E 




(4.64) 



is a correlation matrix. 

a. Theorem 4.5: Schur Recursions 

The quantities a^(k-i-l) and ^(k-N) defined in equations (4.62) and 
(4.63) can be updated according to 



ft N-i(k *0 

4 Li(k-N) 



&(P N + l) 



Q N (k-r 1 ) 

/?N(k -N) 



(4.65) 



and the partial correlation coefficient. /> N - +1 , can be calculated from 

k a S k N (k+ 1) 

" 3s N ( k — N ) 



(4.66) 



b. Proof 

The proof of equations (4.65) follows directly from equations (4.60) 
and the above definitions for ax(k- l) and # N A (k-\). 

Beginning with the definition of the partial correlation coefficient 
given in (4.40). the relationship given in (4.66) is verified as follows 



u 



k t - 1 — ,\ -N 

P N-ri - kv<e kTl r k _^ 



(4.67a) 




b> N (k — N) 



(4.67b) 




b k N : N (k N) 



(4.67c) 
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= a N k - N (k-M)b k %(k-N) 



(4.67d) 



But 



bh-N(k N) - | | y k-N I | [R(k— N)(k-N)jl/2 (4.68) 

and 

/? N k " N (k-N) = [R(x-N)(k-N)]i/2 (4.69) 



Therefore, 



^~ K (k+l) 

fl$“ N (k-N) 



(4.70) 



Using the initial conditions given by (4.59) the following initial 
values are determined for the Schur algorithm 



Q ri*(k~rl) = , T - R (k+1)A 

I I y I I 



a 0 A (k— l)j = r— !R( k+1 )o rO> + Di ... R< k " *) K 

I I y I I 



/?o A (k) = \ R kA 

i y 



j^ 0 A (k) = R k0 R kl • R kK 

: v k 



(4.71a) 



(4.71b) 



(4.72a) 



(4.72b) 



where the parenthesis used to surround the first indices in equation (4.71) simply 
indicate that they are fixed at the value indicated. 

The Schur algorithm implies a filter structure identical to that of 
Figure 4.4. In this case the input vectors are the rows of the correlation matrix 
(normalized by the square root of the diagonal elements). 

5. Synthesis Model 

The original data. V. can be synthesised from the model parameters 
obtained from the Levinson algorithm, equation (4.54). It is also possible to 
regenerate the y A directly from the lattice parameters. The desired result is 
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obtained by solving equations (4.60) for e k * , and r^ 1 . 



-N _ 7~T i2-N + l , k — N 

e k + l - V 1 ~ (PN + i) e k~l PN~l r k-N 



(4.73a) 



r v-V - \/l “ (pn+i) 2 r k-N ~ PN+i e k+V (4.73b) 

Equations (4.73) constitute the synthesis model. They imply a structure 
similar to the analysis model of Figure 4.4, but with the direction of flow of the 
forward error signals reversed. The processing performed at each stage can be 
visualized as in Figure 4.5. A complete third order synthesis model is shown in 
Figure 4.6. Each horizontal path in Figure 4.6 represents a separate synthesis 
model, synthesising a different component of Y. The coordinate system for each 
of these models depends only on values of y A which appear farther down, that is 
they have a smaller value of the index, A. 

Compare Figure 4.3b and 4.5b. It is apparent that the behaviour of the 
backwards error signals is identical in the two cases. Hence, it is possible to 
construct a synthesis model that only reverses the direction of the forward error 
corresponding to the point being predicted. This assumes knowledge of the other 
inputs (zero order forward errors) to the lattice. Such a configuration is 
illustrated in Figure 4.7. 

The amount of knowledge possessed about the signals used for the 
predictions dictates which form of the synthesis model should be used. If little is 
known, then estimates must first be generated which can then be used in the 
prediction. This corresponds to the model of Figure 4.6. If complete knowledge 
is available (either from initial conditions or previous predictions) then Figure 4.7 
can be used. It is also possible to construct models which exploit partial 
knowledge of the input signals and thus fall between these two extremes. In this 
case the known signals should be input as zero order forward errors while t lie 
unknown ones must be estimated. 
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6. Stochastic Fourier Series Interpretation 

From equations (4.72) and (4.40) we can deduce the following result 




(4.74a) 




(4.74b) 



X- 1 



where e k ° +1 is equivalent to y k+1 . These expressions offer an alternate interpretation 
of the lattice filter. Equations (4.74) describe a stochastic Fourier series 
expansion of the forward error sequence where the basis functions are given by 
the backwards error signals. The Fourier coefficients are related to the partial 
correlation coefficients. 

This concludes the review of existing lattice formulations. In the next 
chapter new. multidimensional extensions of this theory are presented. In 
Chapter VI these results are used to derive original nonlinear lattice structures. 
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Figure 4.6: 3-rd Order Generalized Lattice Synthesis Filter 
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V. TWO DIMENSIONAL LATTICE STRUCTURES 



In this chapter a new two-dimensional lattice structure will be derived and 
discussed. Lattice modelling of two dimensional fields has recently received 
considerable attention [Ref. 43,44,45]. In one dimensional lattice modelling, 
updating the order introduces only one new point into the model support. An 
order update in a 2-D lattice model must introduce O(N) new points, where N is 
the model order. Several different solutions have been suggested to this problem. 
The first due to Marzetta [Ref. 43,44], uses a particular ordering of the data to 
reduce the problem to one dimension. He proposed a half-plane support which is 
infinite in one of the two dimensions. This approach, while maintaining several 
of the nice characteristics of 1-D lattices, such as correlation matching and 
producing a minimum phase filter, leads to very long delay filters. 

A different approach, proposed by Parker and Kayran [Ref. 45]. 
simultaneously introduces many points into the support when the model order is 
increased. Their filter uses a quarter plane support and introduces three 
parameters at each order update. Therefore, it lacks sufficient parameters to 
represent all classes of 2-D autoregressive quarter plane filters. More 
importantly, it lacks the property of orthogonality so that the cascading of stages 
does not lead to an optimum filter (better filters are possible using an equivalent 
number of parameters). It's simplicity is attractive and good results have been 
reported using this approach [Ref. 46], 

The theory presented here maintains features of both previous approaches. It 
utilizes the generalized lattice theory presented in Chapter 4 to decompo-' the 
global O(N) point update into O(V) single point local updates. It maintains the 
important property of orthogonality so that the solution at all stages is optimum. 
Although only the quarter plane support case is presented here in detail, the 
theory can be used for any shaped support. It is shown that the Levinson and 
Schur algorithms (see Chapter 4) can be used to solve the 2-D linear prediction 
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problem. In its most general form the lattice contains 0(\ 4 ) parameters while 
there are only 0(N 2 ) points in the filter support. Several structures are presented 
which take advantage of shift-invariance and reduce this requirement to 0(\ s ). 



A. GENERAL FORM OF 2-D LATTICE FILTER 

The theory used is exactly that presented in the previous chapter. The 2-D 
structure results from a careful selection of input data. To illustrate the 
proposed 2-D lattice structure we will consider a 2-D linear prediction problem 
which utilizes a quarter plane support. The 2-D data field is given by 



Y = 




y( A i. A 2)i 



(5.1) 



where (A,,A 2 ) € 2 L K = L K x L K 
where L K is an index set given by 



L K = {0 K 



(5.2) 



Points will be used from this data field in a particular, convenient order. We 
define an ordered index set 



oL h = i(0.0),(l |0),(0,1),(2,0),(0,2), .... (N.O).(O.N). 



(1 : 1), (2.1), (1.2). .... (X.l).(l.N) 

N 

( K — 2 ,K) , ( K l.K l).(K.K- l),(K- l.K),(K.K)j (5.3) 

Other orderings are possible and equally valid. This one is chosen merely to 
illustrate the concepts. The desired ordering of 2 L K . to obtain < 3 L K . can be 
accomplished by the following, computationally efficient algorithm 
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k = 0 

for mO = 0 to K 

for nO = 0 to 2*m0 
if (mod(n0.2) = 0) 
then begin 
i = mO 
j = nO/2 
end 

else begin 
i = nO/2 
j = mO 
end 

0 ’L“ (k) = V (i j) 
k - k + 1 
next nO 
next mO 

In this algorithm qL k (k) has been used to describe the k-th element of the index 
set qL k . while 2 L K (ij) has been used to indicate the (ij)-th element of 2 L K . The 
(K-*-l) 2 elements of 2 L K have been ordered into a one dimensional index set. qL k . 
The elements of qL k can be numbered, consecutively, from 0 to (K-l) 2 -l. The 
notation (k.l) - q will be used to mean: the element of the index set corresponding 
to the q-th element prior to the element (k.l). For example, (2,0)-3 would 
indicate the element (1,0) (see equation (5.3)). Occasionally this notation will be 
abbreviated to simply, kl-q. 

Define the (q-l) order, normalized, forward error associated with the 
prediction of the point y(k.l) from the previous (in the sense of qL k ) (q-l) points, 
as 

euf 1 = a A <V 2 (k.I)y Vj (5.4) 

where the implied summation over (A,.A 2 ) e 2 L K . can be carried out in any order, 
as long as all components are considered. It is preferred to think of (A,. A.,) 
belonging to 2 L K rather than <?L K as this maintains the two-dimensional character 
of the problem. 

The a A q iA '(k,l) can be interpreted as the components of a second order 
covariant tensor. These components, for a range of indices (A,.A 2 ) f (k.l) q (in the 
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sense of (5.3)), or for indices (Aj.A 2 ) > (k.l). are equal to zero. When (A,.A 2 ) = (k.l). 
the component 



&krl(k,1) = I e k r> ! (5 - 5 ) 

where e^ -1 is an unnormalized version of e^ -1 . 

A normalized backward error associated with the prediction of y((k,l)-q) from 
the next q-1 points of qL k , is defined by 

= b^;((k,l)-q)y AlA2 (5.6) 



As with the forward error prediction coefficients, the b A q iA ] ,((k.l) -q) can be 
interpreted as components of a second order covariant tensor. The components 
b> q Ajdk.lj-q) equal zero for the range of indices ( A x , A 2 ) $ ((k,l)-q) or (A,.A 2 ) > (k,l). 
For the case when the index ( A 1; A 2 ) = ((k,l) — q) the component 

b„i:;((k.i)-q) = -rr— (5.7) 

r kl- q 



where is an unnormalized version of r k ^Z‘. 

1. Normalized 2-D Levinson Algorithm 
a. Theorem 5.1 

The 2-D prediction error coefficients can be updated in order using 
the following recursions 



a^ 2 (k.l) 

bqAjfik.l) q) 

where 



©«) 



a TjT J (k.l) 

b^'dfb-D-q) 



@(pq k ‘l 



1 

n / 1 “ [PqY 




and 



p 




l-q- l\ 
r kl— q ( 



(5.8) 



(5.9) 



(5.10) 
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b. Proof (Outline) 

The proof follows directly if the index (k.l) is replaced by a single 
index which runs from 0 to (K-l) 2 -l and thus indexes the elements of qL k . The 
equations (5.8) are identical in form to equations (4.54) and the proof presented 
there can be applied. This approach is equivalent to reordering the data field 
into a vector in the order specified by qL k . 

An alternate proof is possible by generalizing the methods of Chapter 
4. Alternate coordinate systems could be introduced (similar to (4.35) and 
(4.37)). The necessary transformations can be found and all the steps of the 
proof of Theorem 4.2 can be generalized to these higher order objects. These 
concepts will not be explored further here except to note that they could be 
extented to solve problems in any number of dimensions, they are not restricted 
to the two dimensional case being studied here. 

2. Normalized 2-D Error Order Updates 
a. Theorem 5.2 

The two-dimensional prediction errors can be updated in order 
according to the following equation 



— q 




— q— 1 


e kl 


II 

0 
° £ 


e k) 


-q 


— q- 1 


Tkl-q 




r kl — q 



(5.11) 



b. Proof (Outline) 

The proof follows from equations (5.8). If an inner product is formed 
on both sides of the equation (5.8a) with the data field Y equation (5.11a) 
results. Similar arguments can be employed in the verification of (5.11b). 

3. 2-D Form of Schur Recursion 

Define the quantities 



a 



q 1 ’ ] * ( h . 1 } 




(5.12a) 



/C‘h(k,l) — q) 




b^'jid-qJR 



(5.12b) 
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where 



R 





(5.13) 



We note that for the 2-D case the correlation is a fourth order tensor. 

The Schur recursions for the 2-D case are then given by the following 
theorem. 

a. Theorem 5.3 

The generalized 2-D Schur recursions are given by 



A, A 




A,A 0 


a , 1 J (k.l) 


= ©K kl ) 


a q-l (kJ) 


A ,A„ 


A A 


[/V 2 (M-q) j 


£ q -V((k,l)-q) 



and /? q kl can be calculated from 



(5.15) 



b. Proof (Outline) 

The proof follows identical arguments to those of Theorem 4.j and so 
is not given here. 

4. 2-D Lattice Structures 

The derivations presented do not assume shift- invariance. Models are 
built for each point in the data field starting with the point y(K.K) and ending 
with the point y(0.0). All the models are not equivalent, however. In fact, no 
two models are identical. The only model that is quarter plane is the first one. 
that is the model corresponding to the point y(K.K). A quarter plane model for 
any point, y(m.n), in the field can be built by considering an appropriate subset 
of the set Y (equation (5.1)). The subset would ^tart with the point y(m-N.n-N) 
and continue in a quarter plane manner tint il the point y(m.n). for an X-th 
(global) order model. This support can be written 



D 



y ( A i , a 2 



(5.16) 



where 



112 



(A„A a ) € 2 L* n = L m N * L n N 



where 



L 



N 

m 



(m — N),(m 




and 



L 



N 

n 



!(n — N),(n — 




The (N-^l) 2 terms of this index set can be ordered to form 



2y N _ 
0 L mn 



J(m-N,n— N).(mj|N+ l,n — N),(in-N,n— N+ 1), 



(5.17a) 



(5.17b) 



(m- l.n— l),(m,n- l),(m - l,n),(m,n)j (5.18) 

The subset of Y given by (5.16) and the ordering of (5.18) are illustrated in 
Figure 5.1. This could be done for each point in the entire data field Y. yielding 
(K-l) 2 models. If the process is shift-invariant, then all the models would be 
identical. Other simplifications in the model are possible if shift-invariance can 
be assumed. These will be the topic of the next section. In addition, if 
ergodicitv is assumed the required statistical averages can be replaced by 
appropriate spatial ones. 

A second order quarter plane 2-D lattice model for generating the 
prediction error associated with an arbitrary point y(m.n) is illustrated in Figure 
5.2. In this diagram the forward and backward error fields are indicated 
pictorially rather than symbolically. The icons used are defined in the legend on 
the diagram. The large squares, at each stage, indicate the entire support for the 
second order model. The small blank squares indicate the additional support 
(besides the error field squares) used to generate the given error fields. For 
example, the forward error field in the upper right hand square is generated by 
predicting y(m,n) from all the remaining data points in the large square, 
including the one indicated as a backward prediction error. 
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Figure 5.1: Filter Support Ordering (see equation 5.18). 
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The doubly hatched squares, corresponding in this diagram to zero order 
errors, are inputs to the filter. In general, (in later diagrams) doubly hatched 
squares are considered to be inputs, although, they may not always correspond to 
zero order errors. 

The ordering chosen for qL k (equation (5.3)) is only one of many that 
could have been chosen to implement a quarter plane model. This particular 
selection was made for ease of implementation and because it guaranteed that for 
every 2N-fl local updates, a global order update would be completed. In Figure 
5.2 the filter can be visualized as a cascade of increasing order filter sections. For 
every global update. 0(N 3 ) local order updates must be completed. This implies a 
total of 0(N 4 ) updates for an N-th order filter. In general this is too large a 
number to allow these filters to be used for any real time applications. 

In the next section the problem of reducing the complexity of the 2-D 
lattice filter is examined and some solutions are proposed. 

B. REDUCED COMPLEXITY 2-D LATTICE FILTERS 

The assumption of shift-invariance allows certain of the backward prediction 
errors to be considered to be shifted versions of each other. This eliminates some 
calculations. Various structures are possible depending on the types of shifts 
introduced. We note that no single type of shift (neither zf \ z 2 ’L zf'zf ') will 
introduce all the new data elements into the support that are required for a 
global order update. Because of this, additional prediction errors will have to be 
introduced at each stage. This reduces somewhat the advantage gained by the 
shift-invariance assumption. 

Two types of delays will be considered in detail. Initially, several models 
involving diagonal shifts are examined. Later, a model involving a horizontal 
shift is discussed. We begin by introducing a diagonal shift operator. D. This is 
equivalent to multiplication by z, ] z 2 1 in the bivariate /-transform domain. 

Because of the assumed shift invariance the following statement can be 
made 
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R(m.n.m-i.n-j) = Esy(m,n)y(m-i,n-j) < 



(5.19a) 



= EsDjy (m,n) jD[y (m— i,n— j)] 



(5.19b) 



= ESy(m-l,n-l)y(m-i-l,n-j-l) 



(5.19c) 



= R(ij) (5.19d) 

where R() is a correlation function. The correlation is only a function of the 
relative positions of the two points not their absolute positions. If we adopt for a 
moment the Hilbert space formulation (see Appendix A) we conclude that the 
diagonal shift operator is an inner-product preserving operator and so the use of 
the shifted versions of the backward error signals is justified. 

Consider the structure illustrated in Figure 5.3. It represents a third order 
quarter plane lattice filter. At each global stage, (2N-l) 2 -3 lattice coefficients 
are introduced. Therefore, an N-th order model requires 0(N 2 ) parameters. 
Notice that at each stage two new errors are introduced. They each require the 
solution of an (N-l) 2 point prediction problem. For small N this is an 
insignificant number, however, for large N it becomes overwhelming and the 
required number of parameters again becomes 0(N 4 ). It is difficult to analyze this 
structure analytically, as the index sets for each prediction error are different. 
The support for different errors follow different patterns. This, and the 
complexity issue make it a structure that is really only of academic interest. 

The structure of Figure 5.-1 is a true 0(N ? j parameter model. It avoids the 
addition of the new error signals at each stage by introducing them at the outset. 
The structure has a support that differs ^lightly from quarter plane. A more 
significant drawback, however, is that the maximum order of the filter must be 
fixed at the start. If the maximum order is overestimated then some unnecessary 
computations will have been performed. If on the other hand, the maximum 
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Figure 5.3: Reduced Complexity, 2-D Quarter Plane Lattice Filter. 
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given 



required order was underestimated, then a great price must be paid to increase it. 
However, this is considered to be a superior structure because of the regularity 
and complexity reduction it offers. Therefore, a more detailed analysis of this 
model will be performed. 

The ordered index set (equivalent to (5.18)) in this case (for N=2) is given 

by 

oLmn = |(m — 2,n— 2),(m — l,n — 3),(m— S,n — l),(m — l,n — 2),(m — 2,n — 1), 



(m — l.n-l),(m.n-2),(m-2,n).(m,n-l),(m-l,n),(m.n)| 



(5.20) 



The following relations can then be used to advantage to update the backwards 
errors 



**(m— l,n- 1) — ^!^*(m,n)j 


(5.21a) 


r (m-2,n-l) ” ^! r (m-l,n)] 


(5.21b) 


r (m- l,n-2) — ^ r (m,n-l)l 


(5.21c) 


9 9 

r (m-3,n-l) ' ^ r (m-2.n)‘ 


(5.21d) 


**(m- 1 .n 3) — ^(m.n-2); 


(5.21e) 


— 5 ^ _5 

^(m— 2.n-2) ~ .^(m— 1 .n- 1 ) 


(5.21f) 



In general, whenever (m-i.n-j) equals (m.n)-q. then 



(5.22) 



r (m - i- l.i;- j 1 ) r (m -i.n-j) 

Equation (5.22) is a simple rule for exploiting the shift- invariance of the data 
field. 

One final reduced complexity lattice model will be introduced. It will serve 
to illustrate the variety of structures possible and in particular will yield a model 
for which an especially convenient synthesis model can be constructed. The 
model incorporates a horizontal (z,' 1 ) shift, rather than the diagonal shift used in 
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Figure 5.4: Reduced Complexity. 0(N 3 ), 2-D Quarter Plane Lattice Filter 
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the previous models. A horizontal shift also is a inner product preserving 
mapping so that its use is permitted. 

The ordered index set used for this model is given by (for the second order 
case) 

o 2 Lmn = {(m— 2,n— 2),(m— l,n-2),(m,n — 2),(m— 2,n- 1), 



(m-l,n-l),(m,n-l),(m-2,n),(m-l,n),(m,n)| 


(5.23) 


Define a horizontal shift operator, H, by the relation 




y(m.n— l) = H'y(m.n)! 


(5.24) 


The following relations hold due the assumed shift- invariance 
field 


of the given data 


f ( Vi) = H[r,° min) ] 


(5.25a) 


-1 _ IJ-1 

r (m-l,n-l) ^ r (m- l,n). 


(5.25b) 


^(m-2,n-l) — ^ r (m-2.n). 


(5.25c) 


r (m-3,n— l) ^ r (m- 3,n) 


(5.25d) 


— 4 _ II-4 

r (m - 4,n - 1 ) n , r (m- 4,n). 


(5.25e) 


**(m 5.n-l) — ^^(m-5,ri) 


(5.25f) 


In general, whenever, (m-i,n-j) equals (m,n)-q then 




r |m - i.n - j- 1 ) j, 


(5.26) 



T’sing these simplifying relations and equations (5.11). the model of Figure 
5.5 can be deduced. This model still contains 0(N 3 ) parameters, however, the 
actual number is only about one quarter of that required by the previous 
structure (Figure 5.4.) 
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This algorithm shares with the previous algorithm the shortcoming of having 
to estimate the maximum order of the filter at the outset. Despite this, it is 
believed that this model offers a good compromise (probably the best to date) 
between model accuracy and implementation complexity. In the next section it is 
demonstrated that the synthesis form of this algorithm has some particularly 
desirable properties. In Chapter 7, it will be shown that this algorithm is well 
suited for highly parallel VLSI implementation. 

C. SYNTHESIS MODEL 

The synthesis results of the previous chapter are easily extended to the 2-D 
case being considered. The data fields can be regenerated from a knowledge of 
the lattice coefficients and the forward error fields, it is not necessary to explicitly 
calculate the forward and backward error prediction coefficients. 

The desired result is obtained by solving equation (5.11) for e kl 1 and r kl-q- 
This yields 

e„r 1 = s/l - (/>,")* 513 + < (5.27a) 

r kl- q - n / 1 - (Pq M ) 2 r ki-q - Pq kle kl (5.27b) 

These equations establish the method for regenerating the original data field. 
They describe the processing that must be carried out at each stage of the 
synthesis process. As an example of their application, consider the second order 
synthesis model pictured in Figure 5.6. It is the synthesis counterpart of the 
reduced complexity analysis model of Figure 5.5. In order to regenerate the 
original data field processing should be carried out horizontally by rows. For this 
second order model, two rows and two columns of initial conditions must be 
specified. The required zero order forward error sequences will always be 
available from either the given initial conditions or from previous estimates. 

It will be shown that for VLSI implementations it will be more convenient to 
estimate all the zero order forward error sequences. This necessitates that three 
residuals to be input and that all forward error rhannels be reversed. Such a 
structure is illustrated in Figure 5.7. 
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Figure 5.5: Alternate. Reduced Complexity, 0(N 3 ). Quarter Plane Lattice Filter. 
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Figure 5.6: 2-D Quarter Plane Lattice Synthesis Model 
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Figure 5.7: 2-D Quarter Plane Lattice Synthesis Model 
Involving Three Input Residuals 



A synthesis model corresponding to the analysis model of Figure 5.4 would 
have to process the data along diagonals. The required zero order forward error 
sequences would not be available at each stage in this case and so only a model 
analogous to the one of Figure 5.7 can be used. This would require the provision 
of Five input residual signals. 

D. SYSTOLIC IMPLEMENTATIONS 

In the past, most signal processing algorithms were implemented largely in 
software due to their high complexity. More recently, with the advent of VLSI 
technology, there has been a shift towards specialized hardware implementations. 
These offer increased performance at a low per-unit cost. A particularly 
promising class of implementations. First suggested by Rung and Leiserson ■Ref. 

49] . are the so called systolic arrays. These attempt to partition the required 
computations in time and space over an array of identical processing elements, in 
order to increase throughput. 

Rung [Ref. 50: pp. 869] deFmes a systolic array as a network possessing the 
following features: 

a) Svnchrony: The data are rhythmically computed (timed by a global clock) 

and passed through the network. 

b) Regularity (i.e., Modularity and Local Interconnections): The array should 

consist of modular processing units with regular and (spatially) local intercon- 
nections. Moreover, the computing network may be extended indefinitely . 

c) Ternporal Locality: There will be at least one unit-time delay allotted so that 
signal transactions From one node to the next can be completed.* 

dj Pipelinability (i.e.. O(M) Execution-Time Speed-Up): A good measure for the 
efficiency of the array is the following 

~ , T . _ Processing Time in a Single Processor 

Speed t p 1 actor - : — — 

Processing 1 inn* m the Array Processor 

A systolic arrav should exhibit a linear-rate pipelinability . .i.e.. it should achieve 
an O(M) speed-up. in terms of processing rates, when XI h tin number of pro- 
cessor elements (PE’s). 

Methods have been proposed for transforming algorithms into Systolic 
implementations beginning with either an algorithmic description (sec 
Moldovan [Ref. 51]) or a signal- flow-graph (SFG) description (see Rung [Ref. 

50] ). In this chapter we shall be using the second of these methods to transform 



126 



one of the 2-D lattice structures into systolic form. Rather than present a 
detailed discussion of the rules used in this method, we shall simply state them 
and then apply them to produce the desired systolic array. It is hoped that this 
example will clarify the procedures used. If further insight is desired the reader is 
referred to the cited reference. 

1. SFG Transformation Procedure 

The procedure used is based on a cut-set approach. According to Rung 
[Ref. 50. pp. 870] a cut-set is defined as: 

A cut-set in an SFG is a minimal set of edges which partitions the SFG into two 
parts. 

He proposes and proves that the following two rules that can be used to 
transform any computable SFG into a systolic array: 

Rule (i) Time-Scaling: All delays D may be scaled, i.e., D -» D', by a single po- 
sitive integer . Correspondingly, the input and output rates also have to be 
scaled by a factor (with respect to the new time unit D') ... . 

Rule (ii) Delay-Transfer: Given anv cut-set of the SFG, we can group the edges 
of the cut-set into "inbound edges ,f and "outbound edges", depending upon tne 
directions assigned to the edges. Rule (ii) allows advancing k (O') time units on 
all the outbound edges and delaying k time units on the inbound edges, and vice 
versa. It is clear that, for a (time-invariant) SFG, the general system behaviour 
is not affected because the effects of the lags and advances cancel each other in 
the overall timing. Note that the input-input and input-output timing relation- 
ships will also remain exactly the same only if they are located on the same side. 
Otherwise they should be adjusted by a lag of — k’time units or an advance of -k 
time units. 

2. Systolic Implementation of 2-D Lattice Filter 

Using the two rules given in the previous section the 2-D lattice synthesis 
mode! of Figure 5.7 will be mapped into a systolic array. There is some 
flexibility in the design, the result not being unique. The first choice that must 
be made h that of the operation that is to be performed by the basic PE. In the 
case being considered several convenient choices are possible. The simplest 
element would be a multiplier-adder. A slightly higher level operation would be 
that of the single lattice section given by Figure 4.5. A still higher level operation 
is conceivable by grouping several of the lattice sections together. We "hall use 
the second of these choices as it illustrates the general procedure without the 
added complexity inherent in the lower level implementation. 
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The mapping will be done in stages. Initially, the diagram of Figure 5.7 
will be redrawn in an SFG format. Rule (i) will be used to scale all the delays 
appropriately. Then, by successive application of Rule (ii). the SFG will be 
temporally localized (it is already spatially local.) The steps used are outlined as 
follows; 

(1) In Figure 5.8 the algorithm of Figure 5.7 has been redrawn in a SFG form. 
In this and subsequent diagrams delays are indicated by the letter D. The 
lattice section of Figure 4.5, as before, is indicated by the circles at the 
nodes. 

(2) Using rule(i), all delays are scaled by a factor of 6. That is, D -* 6D'. This is 
indicated in Figure 5.9. The input and output signal rates must also be 
scaled by this factor of 6. 

(3) Using the cut-sets indicated in Figure 5.10 the delays are redistributed so 
that temporal locality is achieved. The resulting SFG is indicated in Figure 
5.11. Until now the processing at each node was assumed not to take any 
time. The delays going into each node can be combined with the lattice 
sections and be used to account for the processing time. In this way. the 
structure will appear as in Figure 5.12. In this last Figure, the nodes have 
been shaded to indicate that the operation being performed within them 
consumes one time unit. 

3. Additional Remarks 

In this section we have shown that the 2-D Lattice structures derived in 
this chapter are amenable to a systolic implementation. This is significant as the 
processing of 2-D data Fields such as images in real-time requires high data rates. 
These rates can only be achieved in practice through the use of super-computers 
or specialized hardware'. Due to the high cost of super-computers the second 
alternative is the more practical. With the costs of VLSI production rapidly 
decreasing, it is now cost effective to produce dedicated chips even in very small 
quantities. For large scale productions the cost can be amortized over a large 
number of chips, yielding a low per-unit cost. 
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Although only a single specific implementation has been presented here, 
it indicates the ease with which the other algorithms discussed in this thesis may 
be transformed into forms which can be efficiently implemented in silicon. 

E. SIMULATION RESULTS 

The theory has been proven by computer simulation. Two different order 
models were excited by a white noise process. Several different order estimates of 
each spectrum were generated and compared to the originals. 

1. Example 5.1 

The first model simulated was described by 

y(m,n) - .295y(m-l,n) — ,470y(m,n-l) - 1 - 0.0y(m — l,n— 1) 

— ,055y(m-2,n) -+ .007y(m,n-2) -t- 003y(m — 2,n- 1) 

- .015y(m- 1 ,n— 2) + ,022y(m — 2,n — 2) — u(m,n) (5.28) 

where u(m,n) was a 2-D .zero-mean white noise process. 

The original spectrum is illustrated in Figure 5.13 while the first, second, 
third, and fourth order estimates are given in Figure 5.14. Notice that the 
original model is only second order so that the third and fourth order estimates 
can be used to examine the effects of over modelling. As can be seen from the 
figures the estimated spectra correspond very closely to the originals. Over 
modelling did not noticeably degrade the accuracy of the estimates. 

The actual algorithm used was that of Figure 5.2. Although it is 
unnecessarily complex, it is the most straight forward to implement. The 
generalized 2-D Schur algorithm was used to generate all the required lattice 
parameters. The computer subroutines used to accomplish this simulation are 
included in Appendix B. 

2. Example 5.2 

The second simulation used the following higher order autoregressive 
equation to generate the data 

y(m,n) = — ,47y(m- l.n) — .08y(m.n 1) -+• . 195y(m- l.n — 1 ) 
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Figure 5.8: 2-D Lattice Filter of Figure 5.7 in SFG Form. 
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Figure 5.9: The Delays are Scaled by a Factor of Six. 
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Figure 5.11: Temporally Local Version of 2-D Lattice Filter 
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Figure 5.13: Original Spectrum For Example 5.1. 
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a. First Order Lattice Approximation of Spectrum of Example 3.1. 




b. Second Order Lattice Approximation of Spectrum of Example 5.1. 
Figure 5.14: Lattice Approximations Of Spectrum of Example 5.1 
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c. Third Order Lattice Approximation of Spectrum of Example 5.1. 




d. Fourth Order Lattice Approximation of Spectrum of Example 5.1. 
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— .015y(m,n — 2) -t .055y (m— 1 ,n — 2) - 003v(m- 2,n-2) 



.0067y(m— 3,n) - .015y(m,n — 3) ^ .022y(m-2.n-3) 

-+ .033y(m,n-4) - .085y(m-l,n-4) - .002y(m-4,n — 2) 

~ .0001y(m— 2,n — 4) + .0001y(m — 4,n-4) + u(m,n) (5.29) 

where u(m,n) is a zero-mean white noise process. 

The spectrum of this model is shown in Figure 5.15. The first through 
fourth order estimates of the spectrum are shown in Figures 5.16. The general 
shape of the spectrum is identified in the first order model, although the fine 
detail is not introduced until the fourth order model. The position and relative 
magnitude of the peaks in the spectra are identified with great accuracy in the 
fourth order estimate. 

In the next chapter lattice models are applied to the solution of the 
autoregressive nonlinear modelling problem. 
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Figure 5.15: Original Spectrum For Example 5.2. 
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a. First Order Lattice Approximation of Spectrum oi Example 5.2. 




b. Second Order Lattice Approximation of Spectrum of Example 5.2. 
Figure 5.16: Lattice Approximations Of Spectrum of Example 5.2 
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c. Third Order Lattice Approximation of Spectrum of Example 5.2. 




d. Fourth Order Lattice Approximation of Spectrum of Example 5.2. 
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VI. NONLINEAR LATTICE STRUCTURES 



In the previous two chapters lattice structures have been examined in great 
length and have been applied to the problem of 2-D autoregressive modelling. In 
this chapter we return to the problem of nonlinear autoregressive modelling with 
the hope of solving the problem using a lattice structure. 

There have been several attempts in the literature to fit lattice filters to a 

Volterra like model. The first was due to Parker and Perry [Ref. 15]. They 

proposed a multichannel lattice implementation of an autoregressive-moving 
average nonlinear model. Their model was capable of providing an update in time 
of some terms of the expansion. It, however, does not introduce all the terms that 
are needed to constitute the next higher order model. 

A recent proposal by Zarzycki and Dewilde [Ref. 17] is a true generalization 
of the Levinson algorithm to the case of the autoregressive Volterra type model. 
In their work they considered a non-stationary nonlinear structure and showed 
that the stationary and linear models can be treated as special cases. Their 
model provides a true update in time order, introducing all the required terms. 
They found, as we did with the 2-D model, that not all terms could be 

recursively generated and that some would have to be introduced at each stage 

by other means. For these non-linear models only one error signal must be 
injected at each stage not two as in the 2-D case (if a triangular kernel is used.) 

The model proposed in this chapter is based on the alternate tensor form of 
the nonlinear model developed in Chapter 3. Recall that the model was based on 
interchanging the role> played by time order and nonlinear order. The lattice 
model presented here thus becomes recursive in nonlinear order. This means, for 
example*, that the optimum cubic model is calculated from a knowledge of the 
optimum quadratic model. The theory is based on the generalized lattice 
concepts of Chapter 4 and the 2-D lattice structures developed in Chapter 5. 
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Although only the nonlinear order update model is discussed, some reflection will 
show that a time update can also be performed using the proposed model. 

For simplicity in illustrating the concepts only models involving two delays 
will be considered in this chapter. The theory is equally valid and easily extended 
to more complex situations. 

A. GENERAL NONLINEAR LATTICE MODEL 

As with the 2-D model of the previous chapter, the theory used is that of the 
generalized error order updates given in Chapter 4. A careful choice of the input 
data will yield the desired results. We begin by summarizing briefly the 
nonlinear model to be used, for the case of two delays. 

The estimate of the model output is given by (3.53) 

y(k) - y A '(k 1) • • ■ /"(k-N)H Ai ,. An (6.1) 

To avoid unnecessarily complex algebra we consider the case when N=2. 

y(k) = y Al (k-l)y Aj (k-2)H V2 (6.2) 

The tensor product y Al (k- l)y A2 (k-2) defines a second order tensor, Y(k), with 
components given by 

Y(k) = jy Al (k — l)y A ’(k-2)! = 



1 y (k — 2) y (2) (k-2) 

y(k-l) y (k — l)v(k-2) y(k-1)y' 2 >(k-2) 

y(=)(k- 1 ) y (2) (k 1 )v(k— 2) y (2) (k-l)y |2 )(k- 2) 



y (p) (k-2) 
y(k-l)yW(k-2) 
y (2 >(k-l)y< p) (k — 2) 



y ,pl (k 1) y lp, (k ])y(k 2) y fp) (k - 1 )y ,2, (k 2) 



y (p)(k ])y< p >(k 2) 



(6.3) 

This tensor can be considered to be a shift-varying data field and can thus be 
modelled using the 2-D lattice model of Figure 5.2. A tensor of the form given in 
(6.3) can be formed for each time, k. If the process is time-varying then each of 
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these tensors must be separately modelled. If time-invariance can be assumed, 
then only one model need be developed. In that case the required correlations 
can be calculated over the ensemble of such tensors. 

A shift along the diagonal (or horizontal or vertical ) of the Y(k) tensor is not 
a inner-product preserving operation. This prevents the application of the 
reduced complexity algorithms of Chapter 5 to this nonlinear problem. 

The 2-D lattice model of the tensor Y(k) does not offer a complete solution. 
Notice that y(k), the sample that is to be estimated, does not appear in equation 

(6.3) . Two possible solutions to this problem exist. First, the ordered index set. 

(5.3) , can be extended by one element so that y(k) is included in the model. This 
has the effect of adding a channel to the general 2-D lattice structure of Figure 
5.2. A conceptually different solution is to first model the tensor Y(k) using the 
results of the previous chapter. The backwards error signals that result can be 
used as a basis for a Fourier expansion of the the sample y(k) (see equation 
(4.73) for details.) Both of these approaches lead to identical models but do 
provide alternate interpretations. The second approach will be the one used in 
the derivations of this chapter as it allows the results of the previous chapter to 
be applied with little modification. 

The tensor. Y(k). can be considered to be a two- dimensional data field given 
by 

Y(k) = {y A, (k- l)y >2 (k 2)1 (0.4) 



w here (Aj.A 2 ) L p * L F 

where L p is an index set given by 




(6.5) 



Points will be used from this data field in a particular, convenient order. We 
define an ordered index set 
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(p-pMp-i.pMp.p- 1 )> •••> (o>pMp.°) 



2 

O 



L p = 



(p — l.P— l).(p— 2>p — l).(p— l.p— 2), (O.p-1). (p-1,0), 



...,(1,1),(0,1),(1,0) I (0 > 0)| (6.6) 

The (p+l) 2 elements of 2 L P have been ordered into a one dimensional index 
set, qL p . The elements of qL p can be numbered, consecutively, from 0 to (p+l) 2 -l. 

As in Chapter 5, the notation (m,n) - q will be used to mean: the element of 
the index set corresponding to the q-th element prior to the element (m.n). For 
example, (1.0) - 2 would indicate the element (1,1) (see equation (6.6)). This 
notation will often be abbreviated to simply mn-q. 

Define the (q-1) order, normalized, forward error associated with the 
prediction of the element y m (k— l)y n (k — 2) from the previous (in the sense of qL p ) 
(q-1) points, as 

e^n 1 = a A q ,Aj(m.n)y Al (k— l)y Aj (k — 2) (6.7) 

where the implied summation over- (A,,A 2 ) e 2 L P , can be carried out in any order, 
as long as all components are considered. It is preferred to think of (A,.A 2 ) 
belonging to 2 L P rather than <$1/ as this maintains the multi-dimensional 
character of the problem. 

The '(m.n) can be interpreted as the components of a second order 
covariant tensor. These components, for a range of indices (A,.A 2 ) mn-q (in the 
sense of (6.6)). or for indices (A,,A 2 ) > (m.n) are equal to zero. For the case when 
(A,.A 2 ) = (m.n). the component 

a mn'(m.n) (6.8) 

where e^" 1 is an unnormalized version of e^ n '. 
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A normalized backward error associated with the prediction of v((m,n)-q) 
from the next q-1 points of qL p , is defined by 

fmn-q = b^mn- q)y A ‘(k- 1 )y Aj ( k - 2) (6.9) 



As with the forward error prediction coefficients, the b A q A |(inn-q) can be 
interpreted as components of a second order covariant tensor. The components 
b^Vjmn-q) equal zero for the range of indices (A,,A 2 ) < (mn-q) or (Aj,A 2 ) > (m,n). 
For the case when the index (Aj,A 2 ) = (mn-q) the component 



bmn-q(( m >n)— q) = —r~h (6.10) 

I ' ^mn-q ! ! 

where is an unnormalized version of r q “i q . 

1. Normalized Order Update Recursions 

The following two theorems are stated without proof. The proofs are 
identical to those of Theorems 5.1 and 5.2 of the previous chapter and so are not 
repeated. 

a. Theorem 6.1: Normalized Nonlinear Levinson Algorithm 

The nonlinear prediction error coefficients can be updated in order 
using the following recursions 

a A q ,Aj( m : n ) 
b A q ,A 2 ((m.n)-q) 

where 



( 6 . 



©U' q m 



v ' : (p q mn )-' 



~P q 
1 



= ©(/>c 



a A q Aj( m i n ) 
bA^Ajlmm) 



(6.11) 



‘ r 1 1 \ 



(6.13) 
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b. Theorem 6.2: Normalized Error Order Update Algorithm 

The nonlinear prediction errors can be updated in order according to 
the following equation 



Jnn-q 



= ©(pD 



mn— q 



(6,14) 



In order to introduce the sample y(k) into the model we recognize 
the backwards errors as an alternate coordinate system. In particular the 
following vector is defined. 



jy A : = ! r (o,o )— a 'I 



-0 

r (0,0) 

F ( 0 , 0)-1 
o 

r (0,0)— 2 



F (p-1) 2 -1 
r (0,0)-(p + l) 



2 



-t-1 



(6.15) 



Equations (6.14) correspond to a model structure identical with that 
given in Figure 6.2. In this diagram, the backwards errors given in equations 
(6.15) correspond to those that are shown leaving the uppermost stages of the 
filter. 

The forward error in predicting y(k) from the Y(k) tensor is given by 



e k 



(P- 



= v(k) - Kjry 



.x ' 



( 6 . 16 ) 



Because the backward errors are uncorrelated it is a straight forward 
matter to find the (Fourier) coefficients. K A . They are given by 



K 



f/v (k)rio 



0 0)| 



( b ) r (o.oj— ij 



E 



y(k)r 



(p-i) j -i 

(0.0) (p- 11 = 



t 

'/ 



The m-th component is 



( 6 . 17 ) 
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K m = Esy (k)F(o 0) _„ 



- Ese k m r(o,o)- 



„ 171 r?J~ m ~ m 

e k i M e k r (0,0)-n 



_ m | i k 

ek I Pm 



(6.18a) 

(6.18b) 

(6.18c) 

(6.18(3) 



A normalized version of the forward error defined in (6.16) is thus given by 

m— 1 

E 

A ' = 0 



— m 
e k = 



1 m— i 

|y(k) - E Pvr ( o, 0 )-A-| 1 e k 

! I e k ! I 



(6.19) 



Recognizing that y(k) is a zero order forward error we can write a 
normalized version of the q-th order forward error as 



p q _ 
e k - 






ek 



P* ^o'ol-qj 



( 6 . 20 ) 



This last result follows if (6.19) is iterated, calculating successively 
higher order errors or from equation (4.73) where the equivalence of the lattice 
model and Fourier series was established. 

We conclude that expressing y(k) as a stochastic Fourier series (with 
the backward errors as a basis) amounts to nothing more than the addition of a 
supplementary channel to the existing lattice structure. Equations (6.11) and 
(6.14) hold for this additional channel and can thus be used to update the 
normalized error signals and model parameters. Figure 6.1 illustrates a quadratic 
nonlinear lattice Filter. 

2. Uniqueness Of Lattice Parameters 
a. Theorem 6.3 

The lattice parameterization of the sequence y(k) is unique. That is. 
the lattice parameters given by (6.13) are unique. 
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b. Proof (Outline) 

The truth of this theorem is a consequence of the fact that the lattice 
can be interpreted as a Fourier series expansion of the sequence y(k) in terms of 
the orthonormal backwards errors (see equation (6.17).) The uniqueness of the 
Fourier series is well known [Ref. 47]. 

3. Synthesis Models 

A model analogous to that illustrated in Figure 4.7 can be used to 
regenerate the original sequence. This requires knowledge of the forward error 
residual sequence. The other inputs (zero order forward errors) required can all 
be obtained from past estimates of the sequence or from initial conditions. 

The probability distributions of the output forward error sequences must 
be known if a synthesis model is to be constructed which accurately reproduces 
the original signal statistics. It is not clear that any general statements can be 
made about the nature of these distributions. It has been shown [Ref. 48: p. 357] 
that in the continuous case they will be gaussian and of the same variance as the 
original additive gaussian noise source. In the same reference it is also shown 
that for the discrete case the error sequence will not be gaussian although it will 
be white. The inaccuracies introduced through the use of gaussian noise need to 
be investigated. The stability of these lattice models is also left as an open 
question. 

B. SIMULATION RESULTS 

In order to verify the theory, a computer simulation was performed. Uniform 
white noise was used to excite the system. This choice provides a bounded input 
so that the response of the system used for the simulation could be guaranteed to 
be stable for a suitable choice of system parameters. A sufficient condition for 
the AR model to be stable, given a driving signal whose absolute value is 
bounded on (-a. a), where 0 < a < 1. is given by 

E ' • £ H q a n - a j < 1 (6.21) 

A ,= 1 A n =1 

This condition guarantees that the output never exceeds unity and thus remains 
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Figure 6.1: Quadratic Nonlinear Lattice Filter 
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bounded. This condition is extremely restrictive, however, it does provide a 
simple rule for building models which can be used for simulation purposes. (The 
issue of the inaccuracy introduced by using uniform noise to drive the system has 
been avoided.) 

The model tested was 



[Ha,aJ 



-.1 .22 .02 
.02 .09 .001 
.2 .05 -.03 



(6.22) 



Using the nonlinear form of the Levinson algorithm, equation (6.11) the 
following model parameters were estimated for two repetitions of the experiment, 



-.1003 


.26038 


.04568 


.01402 


.10198 


.03877 


.1998 


.03230 


-.03965 


-.1009 


.2612 


.06213 


.01151 


.09917 


.03125 


.19020 


.04199 


-.01432 



( 6 . 23 ) 



( 6 . 24 ) 
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VII. CONCLUSIONS AND DISCUSSION 



A. SUMMARY OF NEW RESULTS 

In this dissertation the use of tensor concepts in the field of signal processing 
was investigated. The research was successful in a number of areas, extending 
known results and introducing some new ones. In particular, tensors were used 
to study nonlinear and multidimensional systems. 

The nonlinear modelling problem was formulated using tensors. By 
interchanging the roles played by the time order and the nonlinear order a new 
form, different from (although equivalent to) the traditional Volterra Series was 
developed. Using this new form, tensor equivalents of the discrete-time Wiener- 
Hopf and Yule- Walker equations were derived. These equations can be solved for 
the unknown system parameters. Conditions for the solution of the Wiener-Hopf 
equations, without requiring matrix inversion, were specified. This resulted in a 
computational saving of 0((p+l) 5(N "^) operations, where N is the largest time 
delay present and p is the highest exponent present in the system model. It was 
further shown that linear adaptive algorithms, such as LMS and RLS, can be 
applied to solve for the system parameters. 

Existing Lattice filter algorithms were reformulated in a tensor framework. 
It was shown that they can be considered to be equivalent expansions in 
alternate coordinate systems. These results were then applied to the solution of 
the 2-D autoregressive modelling problem. Several new 2-D lattice structures 
were deduced. These models are not efficient in the sense that an AR model 
possessing 0(N 2 ) parameters would require 0(V) parameters when recast into a 
lattice form. It was shown that with proper assumptions of shift-invariance the 
complexity of the lattice models can be reduced to 0(N ? ') parameters. 

The 2-D lattice models developed were then used to deduce a nonlinear 
lattice model. This model differs from that of other researchers in that it allows 
updates to be performed in the nonlinear order as well as time order. 
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Finally, it was shown that these lattice models are amenable to systolic array 
implementations. 

B. FUTURE DIRECTIONS FOR RESEARCH 

The multidimensional and nonlinear lattice theories are by no means 
complete. So far, several new lattice models have been developed. It now 
remains to investigate the properties of the models. 

For the 2-D lattice Filters, conditions for model stability have yet to be 
established. The stability of multidimensional systems is a much more difficult 
problem than for the 1-D case because it is usually impossible to factor 
multivariate polynomials. However, necessary and sufficient conditions for AR 
2-D model stability have been established. It is believed that these can be used 
to derive lattice stability conditions. 

Another important topic is the synthesis of 2-D transfer functions using 
lattice structures. Stated differently, the problem is to design a 2-D lattice 
parameter filter that will implement a given 2-D frequency response 
characteristic. 

For the nonlinear lattice structures similar issues need to be addressed. In 
this case, however, the problems are more complex since the stability and 
synthesis questions have not been solved for the nonlinear AR models. 

Before the synthesis problem for the lattice implementation can be solved 
several more basic questions need to be answered. The first is the definition of a 
useful and meaningful specification of the desired filter characteristic. The 
nonlinear AR filter affects not only the frequency content of the driving signal 
but also its probability distribution. The filter will also have different frequency 
characteristics for different driving signals. All of these effects should be included 
in the specification. 

It has been shown that for the continuous case, a nonlinear whitening filter 
will yield a signal with a gaussian probability distribution [Ref. 48: p. 357 . For 
the discrete case it has been shown that this is not true. The inaccuracy 
introduced by approximating the input to the synthesis model by a gaussian 
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introduced by approximating the input to the synthesis model by a gaussian 
signal is an important question. 

Stability of the nonlinear lattices is also an important question. No simple 
stability tests exist for the nonlinear AR models. Their recursive nature makes 
stability difficult to analize. In general, stability will be a function of the input 
as well as the system which significantly complicates the problem. 
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APPENDIX A: ALTERNATE PROOF OF THEOREM 4.2 



In this appendix an alternate proof of the generalized error order update 
equations (4.40) is provided. The proof presented here relies on geometric 
arguments which are possible if an abstract mathematical framework is adopted. 
In this appendix, a random process will be considered to consist of a time-series 
of random variables. Each of these random variables is thought of as a vector in 
an infinite dimensional Hilbert space. For a rigorous discussion of these concepts 
see for example [Ref. 51,52]. 



A. DEFINITIONS AND FORMULATION 

A discrete random process Y = jy(i): 0 < i < k| can be considered to span a 

K + l dimensional subspace S 1 * 41 of an infinite dimensional Hilbert space S 
(assuming completeness and the linear independence of the y(i)). The inner 
product on this space is defined to be 



<u,v> = 




(A.l) 



where u and v are elements of S. This inner product induces a norm 
u = <u.u>’ 2 



(A.2) 



The error. e/L , . in predicting the element y(k+l) from the previous N 
elements of Y is given by 



<V-i = £ h A N (k- l)y(A) 

A- 0 

where 

h N (k i) o o h k N N ,, h k N N _ 2 h, N io o 

A normalized version of the forward error is given by 



(A.3) 



(A.4) 



(A. 5a) 



1G0 



(A. 5b) 



= £a A N (k-l)y(A) 



A=0 



where 



a, N (k+l) 



h A N (k+l) 



(A.6) 



I I e k+l i I 

The backwards prediction error, r^, is the error associated with the 
prediction of y(i) given the next N elements of Y. It is given by. 



r k N -N= E h* N (k-N)y(A) 

A=0 



(A.7) 



where 



fh/*(k-N); = jo - • - o i -h k N . N+1 -h k N . NT2 • ■ • 
A normalized version of can be defined as 



-K o 



(A.8) 



r k-N 



r N 
r k-N 



r N 
^k- N 



(A. 9a) 



E b A N (k-N)y(A) 

A = 0 



(A. 9b) 



where 



b^(k-N) -- 



hjT(k-N) 

r k-N i I 



(A. 10) 



By orthogonal we mean that given two elements of the space S. u and v. say. 



then 



u. v - - 



0 for u v 

u ! for u = v 



(A. 11) 



This will be indicated symbolically as 



(A. 12) 



The expression 



161 



(A. 13) 



u _ S M 

implies that u is orthogonal to all elements of the subspace S M . 

The symbol ' ©‘ is used to indicate direct sum. For example 

S 3 = S' © S 2 (A. 14) 

means that the space spanned by S 3 is the the space spanned by the Cartesian 
product of the underlying sets of the spaces S 1 and S 2 [Ref. 47: p. 196]. 

Finally, the symbol ’ v ’ will be used to mean the span of. 

These definitions are sufficient to state and prove the generalized error order 
update theorem. 



B. ERROR ORDER UPDATE RECURSIONS 
1. Theorem 4.2 

The (N + l) order errors can be calculated from the N- th order errors 
through the recursion 



-N + l 

e k + l 



— N + l 
r k-N 



1 


1 -Pn+i 






\/l - (Pn*i) 2 


_ Pn + i 1 




-N 

r k-N 



(A. 15) 



where 



P.N + 1 = < e k + 1 - r k^-N> (A. 16) 

is called the partial correlation coefficient or reflection factor [Ref. 17: p. 
37]. 

2. Proof 

We have that 



\ 



v jy(k 



N-l) >(k N-2) 



blit 



i*k- \ 



\ t ] 



> ( k N) >(k N 1) 



Therefore. 



> (k)j' 

y(k)j 



(A. 17) 



(A. 18) 
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S k N "> - v F k N ; N y(k-N-l) • • • y(k) 



(A. 19) 



which can be written as 



S k N+1 = S k N 0 v r k % 



(A. 20) 



For the updated forward error, e k ”V 5 the following is true 



N + l | q N V 1 

e k + l i_ 1 r k-N| 



(A.21) 



Since ejl, _ S k N we need only orthogonalize it (by a Gram-Schmidt procedure) 
with respect to v ^v. N |- to obtain eJVi 1 - Thus. 



„ N + l _ N /^-N 
e k + i - e k -M — t^r k _] 



where 



N — N 

<e kTl' r k- N > 

- N 
r k- N 



(A. 22) 



(A. 23a) 






e k*-l 



(A. 23b) 



Therefore, 

N 

-N-M -N k — N 

e k-^l r ~ e k-l ~ r k \ 

e k-l 



(A. 24) 



Proof of the fact that 

e k N -) l 

’k N - 1 ' \/i (^n*i) j 



(A .25) 



will not be repeated here (see Chapter 4. Section B). 

Similar arguments can be used to verify the backward error order update 
equation given in (A. 15). 
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APPENDIX B: FORTRAN PROGRAM LISTINGS 



**X***X*****xx*****:»:******x**************xx*:i:***x*x**x**********x****** 



2DLAT MAIN PROGRAM 

PROGRAM TO TEST 2D LATTICE SUBROUTINES 
WRITTEN 29 APRIL 1985 

*********************************************************************** 



REAL'S ALPHA(26,26),BETA(26,26),A(26,26),B(26,26),RHO(26,26) 

REAL'S HZ(25),Y(256,256),R(26,26) 

REALM GAIN(40,40) 

DEFINE THE AR FILTER COEFFICIENTS 
DATA HZ/1.0, -.23, .12, .111, 21*0.0/ 

DATA HZ/l.O,-. 470, .007, .295, 0.0, .015, .055, .003, .022, 16*0.0/ 

DATA HZ 1 0, .03, - 015, -.011, .033, -.47..195, .055.0.0, -.085, 0.0. 0.0, 

*.003,. 022.-.0001..0067, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -.002, 0.0, .0001/ 

DEFINE OTHER PARAMETERS 

MODEL SIZE 

N = 5 
MN = N*N 

ACTUAL SYSTEM SIZE 
NA = 5 

MNA = NA*NA 
C 

C NUMBER OF POINTS IN THE PLOT IS IP X IP. 

C 

IP = 40 
C 

C NUMBER OF DATA POINTS TO ESTIMATE CORRELATIONS IS IYS X IYS. 
C 

IYS - 128 
C 

C CALL APPROPRIATE SUBROUTINES 
C 

( A L i. TX FUN ( H Z ,N A . G A I N .IP.0) 

CALL PLOT3((.AIN IP IP ORIGINAL SPECTRUMS) 

C A L L G E N R A T ( H Z . N A . Y , I Y S ) 

CALL CORLAT(Y,IYS,R,N) 

CALL SC'HURfRHO.R. ALPHA. BETA. MN) 

CALL LEYSON) A.B.RHO.R .MN) 
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DO 10 I = l, MN 
HZ(I) = A(l.I)/ A(I,I) 

10 CONTINUE 

CALL TXFCN(HZ,N.G AIN.IP.l) 

CALL PLOTS (GAIN. IP. IP, '4-TH ORDER LATTICE APPROXIMATIONS') 

STOP 

END 
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non 



***:fc*st:3if****:<c**:r*x**:X5|£X*xx**:rx:fcxxj|:sfcxxx>'Sir:*:r*xj(c:fc 



* x* x 



q************* 

c 

C SUBROUTINE GENRAT 
C 

C THIS SUBROUTINE GENERATES A 2D DATA FIELD FROM THE AR 
C FILTER COEFFICIENTS. A WHITE NOISE INPUT AND WHITE NOISE 
C INITIAL CONDITIONS ARE USED. 

C 

C WRITTEN 29 APRIL 1985 

C 

£«*********** ****************************** ********* ********** *********** 
C 

SUBROUTINE GENRAT(HZ,N,Y,IYS) 

REAL*8 HZ (25 ),Y( 256.256), ADD 
C 

C FETCH THE RANDOM NUMBER GENERATOR SEED 
C UNDER FORTHX. STORE SEED IN FILE 13. 

C UNDER DISSPLA ENTER SEED EXPLICITLY. 

C 

C READ(13,10) IY 
CIO FORMAT(IlO) 

C REWIND 13 
IY = 817346 

GENERATE THE RANDOM FIELD 

MN = N*N 
MNM 1 = MN- 1 
DO 40 I = l.IYS 
DO 30 J - l.IYS 

Y(I.J) - URAND(IY) - .5 
KI = 0 
KJ — 0 

DO 20 K =■ l.MNMl 
IF (MOD(K.N) NE.O) GOTO 15 
KI — KI — 1 
KJ 0 
GO TO 16 

15 KJ = KJ - 1 

16 ir ((II-KI) LEO) OR.f(J-KJ).LE.O)) GOTO 17 

ADD X (I-KI.J-KJ) 

GO TO 18 

17 ADD UR AND (I Y) - .5 

18 Y(I.J) Y(I.J) - HZ( K • 1 )’ ADD 

20 CONTINUE 

30 CONTINUE 

40 continue 

c 

(' STORE THE RANDOM NUMBER SEED 
C 

C WRITE(13.50) IY 

C50 FORMAT(IIO) 

C 
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RETURN 

END 
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** * * 



**y»i***!«*l******jf****x*?4****?**» 



****acj|c*******:***±******!*:**x*:*ar:r* 



SUBROUTINE CORLAT 

THIS SUBROUTINE PRODUCES A CORRELATION MATRIX FROM 

A 2-D DATA FIELD IN AN ORDER 

WHICH IS COMPATIBLE WITH SUBROUTINE SCHUR 

WRITTEN 30 APRIL 1985 



SUBROUTINE CORLAT (Y,IYS,R,N) 
REAL*8 Y(256,256),R(26,26),SUM 

DEFINE CONSTANTS 

MN = N*N 
IR = 0 



BEGIN OUTER LOOP 

DO 100 MPl = l.N 
MO = MPl - 1 
LLIM = 2*MPl - 1 
DO 90 L = l.LLIM 
L0 = L - 1 
II = MO 
Jl = LO/2 

IF (MOD(LO,2).EQ.O) GO TO 10 
II = Jl 
Jl = M0 
10 JR = IR 

IR = IR + 1 
KBOT = L 
DO 80 NPl = MPl.N 
NO = NPl - 1 
KLIM 2* NPl - 1 
DO 70 K = KBOT. KLIM 
K0 ^ K-l 
I = NO 
J = K0 2 

IF (MOD(K0.2).EQ.0) GO TO 20 
I J 
J = NO 

20 IOFF = 1-11 

JOFF = J - Jl 
KlMAX = 1YS 
K 1 MIN = IOFF - 1 
IF (IOFF.GT.0) GO TO 30 
KlMAX = IYS t- IOFF 
KlMIN = 1 

30 K2N1AX - IYS 
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K2MIN = JOFF -r 1 
IF (JOFF.GT.O) GO TO 40 
K2MAX = IYS - JOFF 
K2MIN = 1 
40 SUM = 0.0 

JR = JR + 1 

C WRITE(6,44)IR,JR,I1,J1 ) I,J,10FF,J0FF 

C44 F0RMAT(8(2X,I4)) 

DO 60 Kl = K1MIN,K1MAX 
DO 50 K2 = K2MIN,K2MAX 
SUM = SUM + Y(Kl,K2)*Y(Kl-IOFF,K2-JOFF) 

50 CONTINUE 

60 CONTINUE 

R(IR,JR) = SUM/(KlMAX-KlMIN+l)/(K2MAX-K2MIN^l) 
70 CONTINUE 

KBOT = 1 
80 CONTINUE 

90 CONTINUE 
100 CONTINUE 
C 

C FILL IN THE SYMMETRIC HALF OF CORRELATION MATRIX 
C 

DO 120 I = 2.MN 
IMl = 1-1 
DO 110 J = l.IMl 
R(I.J) = R(J.I) 

110 CONTINUE 
120 CONTINUE 
C 
C 

RETURN 

END 
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c 

C SUBROUTINE TXFCN 
C 

C GENERATES A 2-D FREQUENCY RESPONSE GIVEN THE TRANSFER 
C FUNCTION COEFFICIENTS 
C 

C WRITTEN BY DR. B. MADAN 
C MODIFIED BY P.J. LENK 29 APRIL 1985 
C 

C PARAMETER KO INDICATES THE ORDERING OF THE COEFFICIENTS 
C - KO = 0: PARAMETERS IN ROW-MAJOR ORDER 

C - KO = 1: PARAMETERS IN TWIDDLED ORDER 

C 

q****************************** ****************************** *********** 

c 

SUBROUTINE TXFCN(HZ,N,GAIN,IP,KO) 

REAL*8 HZ(25),CONVRT(5,5) 

REALM G AIN(40,40) 

COMPLEX CSUM 
COMPLEX *8 ARG1.ARG2 
C 

C DEFINE CONSTANTS 
C 

PI - 3.14159 

DETERMINE ORDER AND REORDER IF NECESSARY THE COEFFICIENTS 

IF (K0.EQ.0) GO TO 60 
JR - 0 

DO 30 MPl - l.N 
MO = MPl - 1 
LLIM = 2*MPl - 1 
DO 20 L = l.LLIM 
L0 = L - 1 
I = MO 
J = L0 2 

IF (MOD(L0.2).EQ.O) GO TO 10 
I = J 
J = M0 
10 JR JR - 1 

CONVRT(l- l.J-1) - HZ(JR) 

20 CONTINUE 
30 CONT1N1 L 
( 

C TRANSFER THE COEFFICIENTS BACK TO HZ 
C 

JR - 0 

DO 50 I = l.N 
DO 40 J = l.N 
JR = JR » 1 
HZ(JR) = CONVRT(I.J) 

WRITE(6.35)HZ(JR) 
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35 FORM AT (2X, El 2.5) 

40 CONTINUE 

50 CONTINUE 
C 

C PROCEED WITH TRANSFER FUNCTION EVALUATION 
C 

60 DW = 2.0*PI/FLOAT(IP-1) 

DO 100 IWl = i,IP 
W1 = DW * (IWl - 1)- PI 
DO 90 JW1 = 1,IP 
A1 = Wl 

W2 = DW * (JWl - 1)- PI 
IZ = 0 

CSUM - CMPLX(0. 0,0.0) 

DO 80 I = 1,N 
ARGl = CMPLX(0.0,-A1) 

A1 = A1 + Wl 
A2 = W2 
DO 70 J = l.N 
IZ = IZ + 1 

ARG2 = CMPLX(0.0,-A2) 

A2 = A2 + W2 

C WR1TE(6, 77)1, J, ARGl, ARG2,IZ,HZ(IZ), CSUM 

C77 FORMAT(2(2X.l3),4(2X,El2.5),2X,l3,3(2X,El2.5)) 

CSUM = CSUM + CMPLX(SNGL(HZ(IZ)),0.0)*CEXP(ARG1+ARG2) 
70 * CONTINUE 

80 CONTINUE 

GAlN(lWl.JWl) - 1.0/CABS(CSUM) 

GAIN(IWl.JWl) = GAL\(IWl,JWl)*GAIN(IWl,JW 1) 

C WRITE(6.78)IWl.J Wl.CSUM.G AIN(IWl.JWl) 

C78 F0RMAT(2(2X.I3).3(2X E12.5)) 

90 CONTINUE 

C WRlTE(6,]8)IWl.(GAIN(IWl,I),I=l,IP) 

C18 FORMAT( 1X.I3,5(2X.E12.5)) 

100 CONTINUE 
RETURN 
END 
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********* 



* ** * * * 



* * 



************************** 



******************** 



SUBROUTINE PLOTS 

THIS IS A ROUTINE TO USE THE DISSPLA PACKAGE TO 
DRAW A THREE DIMENSIONAL PLOT OF A 2-D FILTER’S 
FREQUENCY RESPONSE. 

WRITTEN BY DR. B. MADAN 
MODIFIED BY P.J. LENK 29 APRIL 1985 



*********************************************************************** 



SUBROUTINE PLOT3(A,IM,JN, LABEL) 
DIMENSION A(IM,JN) 

INTEGER LABEL) 100) 

CALL TEK618 

INITIALIZE THE PLOTTING DEVICE 

WR1TE(6,51) 

51 FORMAT)’ 1 CALL TEK618 OK’) 

C CALL THE DEVICE 

CALL RESET(’ALL’) 

WRITE(6,52) 

‘ 52 FORMAT)' 2 RESET ALL OK’) 

C CALL SETUP FOR CUBE 

A 1 = FLOAT(IM) 

A2 = FLOAT(JN) 

CALL PAGE) 1 1 .9.5) 

WR1TE(6,53) 

53 FORMAT) ' 3 CALL PAGE OK’) 

C CALL PAGE(A1,A2) 

C CALL PHYSOR(0.,0.) 

CALL AREA2D(7. 0.7.0) 

WR1TE(6.54) 

54 FORMAT)' 4 CALL AREA 2-D OK’) 

CALL S1MPLX 

CALL HEIGHT). 2) 

CALL HEADINfL ABEL. 100.4.1) 

CALL HEIGHTIO. 14) 

CALL VIEW (-10.0.-5.0.15.0) 

CALL Vf »LM3D( 12 . 0.1 2.0. 12.0) 

( C\L1. X \XI'' LABELLING ROUTINE 

CALL X3NAMI.C \\ 2 f.3) 

C CALL 'i AXIS LABELLING ROUTINE 

CALL YSNA.MLC Wl 8'.3) 

C CALL Z AXIS LABELLING ROUTINE 

CALL Z3.NAME) ".2) 

C CALL THE SURFACE PLOT ROUTINE 

C 

CALL GRAF3D(-E0.0.2,E0.-E0.0. 2. 1.0.0.0.2.0.16.) 
1500 CALL SVRMAT) A . I .IM I.JN.O) 
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CALL ENDPL (1) 
CALL DONEPL 
RETURN 
END 
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SUBROUTINE SCHUR 

CALCULATES THE REFLECTION FACTORS FROM THE CORRELATION MATRIX 
WRITTEN 29 APRIL 1985 

*********************************************************************** 



SUBROUTINE SCHUR(RHO,R, ALPHA, BETA, N) 

REAL*8 RHO(26,26),R(26,26),ALPHA(26,26),BETA(26,26),RNORM,T 

INITIALIZE THE ALPHA AND BETA ARRAYS 

DO 10 I = 1,N 
DO 5 J = 1,N 

ALPHA(I.J) = R(I,J)/DSQRT(R(I,I)) 

BETA(l.J) = ALPHA(LJ) 

RHO(I.J) = 0.0 
CONTINUE 

WR1TE( 16,7) (ALPHA (I. J),J=1,N) 

FORMAT(5(2X,El2.5)) 

CONTINUE 

BEGIN CALCULATING THE REFLECTION FACTORS 

DO 50 J = 2.N 
NJl = N - J -i- 1 

DO 40 1 = l.NJl 
Jll = J — 1 - i 
1P1 = 1 - 1 

RHO(l.Jll) = ALPHA (I,JIl) /BETA (1P1.JU) 

RNORM = DSQRT( 1.0 - RHO( I . JI 1 ) *RHO ( I. J 1 1 )) 

DO 50 K = l.N 
T = ALPHA(l.K) 

ALPHA(l.K) - ( ALPHA (l.K)-RHO(l. Jll) 'BETA (IPl.K)), RNORM 
BETA(l.K) - (BETA(lPl.K)-RHO(l.Jll)*T) RNORM 
50 CONTINUE 

40 CONTINUE 

C \\ R ITE( 16.42 )J.(( ALPHA (I.K).K 1 ,N).l l.N) 

(.'•42 FORMAT) 2X.I5.4(2X.E12.5)) 

C WRITE) I6.42)J.((BETA(I.K).K - l.N). I- l.N) 

50 CONTINUE 

C 

C 

RETURN 

END 
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Q*)ie^)(tX!|C******|!*!(!i»¥***********’ r *^*******)l t » * * * * * » * 

C 

C SUBROUTINE LEVINSON 
C 

C GENERATES THE AUTOREGRESSIVE MODEL COEFFICIENTS GIVEN THE 
C REFLECTION FACTORS 
C 

C WRITTEN 29 APRIL 1985 
C 

Q*** ******** ************ ********* *************************************** 

c 

SUBROUTINE LEVSON (A,B,RHO,R,N) 

REAL*8 A(26,26),B(26,26),RHO(26,26),R(26,26),RNORM,T 
C 

C INITIALIZE THE A AND B MATRICES 
C 

DO 20 I = 1,N 
DO 1 0 J — 1,N 
A (I, J ) = 0.0 
B ( I, J) = 0.0 
10 CONTINUE 

A (1,1) = 1.0/DSQRT(R(I,I)) 

B(I.I) = A (1,1 ) 

20 CONTINUE 
C 

C CALCULATE THE AR PARAMETERS USING LEVINSON’S ALGORITHM 
C 

DO 60 J = 2,N 
NJl = N - J + 1 
DO 50 I = I, NJl 
1J1 = I - J - 1 

RNORM = DSQRTfl.O- RHO(LIJl)*RHO(I.IJl)) 

DO 40 K = FIJI 
T = A(I.K) 

A(I.K) = ( A(I.K) - RHO(I.IJl)*B(I+l.K)) 'RNORM 
B(I.K) = (B(l-l.K) - RHO(I.IJl)*T) 'RNORM 
40 CONTINUE 

50 CONTINUE 

C WRITE) 16.77) J.((A(1. K).K = 1 ,N).I= l.N) 

C WRITE(16.77)J.((B(I.K).K = l.N). I -l.N) 

C'77 FORMAT!. 2X.I3.4 (2X.EI 2.5) ) 

60 CONTINUE 
DO 66 I = l.N 
RNORM = A(l.I) A(l.l) 

C WRITE((>.65)RNORM 

c WRITE) 1 7, 65 )R NORM 

( '65 FORMAT(2X.El2.5) 

66 CONTINUE 

C 

C 

RETURN 

END 
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* * * ****************** ± X ******** ± x x ***** X * * * * X X * X * * * X x X X * X * *■ X x x ± * ** yp. 

c 

C FUNCTION URAND 
C 

C TAKEN FROM "COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS" BY 
C G.E. FORSYTHE. M.A. MALCOLM. AND C.B. MOLER 
C PRENTICE-HALL, ENGLEWOOD CLIFFS, NJ., 1977, P. 246. 

C 

C 

REAL FUNCTION URAND (IY) 

INTEGER I A, IC, IT WO, M2, M, MIC 

C URAND IS A UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND 
C SUGGESTIONS GIVEN IN D.E. KNUTH (1969), VOL. 2. THE INTEGER IY 
C SHOULD BE INITIALIZED TO AN ARBITRARY INTEGER PRIOR TO THE FIRST 
C CALL TO URAND. THE CALLING PROGRAM SHOULD NOT ALTER THE VALUE OF IY 
C BETWEEN SUBSEQUENT CALLS TO URAND. VALUES OF URAND WILL BE RETURNED 
C IN THE INTERVAL (0,1). 

C 

DOUBLE PRECISION HALFM 
REAL S 

DOUBLE PRECISION DATAN,DSQRT 
DATA M2/0/.ITWO/2/ 

IF(M2.NE.O) GO TO 20 

IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH 



M=1 

10 M2=M 
M=ITWO*M2 
IF(M.GT.M2) GO TO 10 
HALFM=M2 



COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD 

IA = 8 x IDINT(HALFM*DATAN(l.D0)/8.D0)-5 
IC — 2*IDINT(HALFM*(.5D0-DSQRT(3.D0) 6.D0))+ 1 
MIC (M2-IC)-M2 

C 

C S IS THE SCALE I ACTOR FOR CONVERTING TO FLOATING POINT 
C 

S=.f» HALFM 
C 

C COMPETE NEXT RANDOM NTMBER 
(' 

20 n IN’ * I A 

c 

(■ THE FOLLOWING STATEMENT IS FOR COMPUTERS W HICH DO NOT ALLOW 
C INTEGER OVERFLOW ON ADDITION 
C 

IF(IY.GT.MIC) IY=(IY-M2)-.\12 
C 

IY IN' — K ’ 
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THE FOLLOWING IS FOR COMPUTERS FOR WHICH THE WORD LENGTH 
FOR ADDITION IS GREATER THAN FOR MULTIPLICATION 

IF(IY/2.GT.M2)IY=(JY-M2)-M2 

THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER OVERFLOW 
AFFECTS THE SIGN BIT 

IF(IY.LT.0)IY=(IY+M2)+M2 

URAND=FLOAT(IY)*S 

RETURN 

END 
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c 

C PROGRAM NLMAIN 
C 

C THIS IS THE PROGRAM TO TEST THE NONLINEAR LATTICE MODEL 

C 

C WRITTEN 7 MAY 1985 
C 

Q* ************************************* *************** ****************** 

c 

REAL*8 ALPHA(26,26),BETA(26,26),A(26,26),B(26,26),RHO(26,26) 

REAL*8 HZ (5, 5), Y( 10000), R (26, 26) 

DEFINE SYSTEM PARAMETERS 

DATA HZ/. 2, .02, -.1.0. 0,0.0, .05, 0.09, .22, 0.0, 0.0, -.03, .001, 0.02, 0.0 

*, 0 . 0 , 10 * 0 . 0 / 

DEFINE SYSTEM CONSTANTS 

DO 5 I = 1,5 
K = 6 - I 

\VR1TE(6,4)(HZ(K,J),J=1,5) 

FORMAT(5(2X,El2.5)) 

CONTINUE 
IYS = 1000 
NA = 3 

MNA - NA * NA 
MNAPl = MNA + 1 



DEFINE MODEL CONSTANTS 
N = 3 

MN = N * N 
MNPl = MN - 1 
C 

C CALL SUBROUTINES 
C 

CALL NLGENf Y.N.HZ.IYS) 

CALL NLCLAT(Y.IYS.R.N) 

DO 30 I 1 .MNPl 

W R ITE( 16.20) (R(I.J). .1 — 1, MNPl) 

20 FORMAT) ,5(2X.El2.5)) 

50 CONTINUE 

CALL SCHUR(RHO.R. ALPHA. BETA. MNPl) 
DO 60 I 1 .MNPl 
\\ R ITE( 16.20) (R HO (I. .]),.! — 1 .MNPl ) 

60 CONTINUE 

CALL LEVSON) A,B,RHO,R ,.\lNPl ) 

DO 70 I = 1 .MNPl 
WRITE) 16, 20)(A(LJ),J=1, MNPl) 

70 CONTINUE 
C 



178 



STOP 

END 
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C 

C SUBROUTINE NLGEN 
C 

C THIS SUBROUTINE GENERATES AN OUTPUT SEQUENCE FROM THE SYSTEM 
C DESCRIBED BY THE MODEL PARAMETERS CONTAINED IN H(,). IT 
C USES WHITE NOISE UNIFORM ON (-.5, .5) TO EXCITE THE SYSTEM. 

C THE INITIAL CONDITIONS ARE ALSO DRAWN FROM THIS DISTRIBUTION. 

C 

C WRITTEN MAY 7 1985 
C 

C 

SUBROUTINE NLGEN(Y,N,H,IYS) 

REAL*8 Y(10000),H(5,5),DSEED 
C 

C FETCH THE RANDOM SEED 
C 

READ(13,10) IY 
10 FORMAT(IlO) 

REWIND 13 

C DSEED = DFLOAT(IY) 

C 

C SET UP THE INITIAL CONDITIONS 
C 

Y ( 1 ) - 2.*(URAND(IY) - .5) 

Y(2) = 2.*(URAND(IY) - .5) 

C CALL GGNML(DSEED,IYS,Y) 

C 

C CALCULATE THE REMAINING VALUES OF THE SEQUENCE 
C 

DO 40 I = 3.IYS 

C Y(I) = 2.*(URAND(IY) - .5) 

C Y( I) = ,2*Y(1) 

Y(I) = URAND(IY) 

C DO 30 J = l.N 
C JM1 = J - 1 

C DO 20 K l.N 

C KMl = K - 1 

C Y(l) Y(I) - H(J,K)*COORD(Y(I-l ).JMl)”COORD( Y(l-2).KMl ) 

( '20 CONTINUE 

C30 CONTINUE 
10 CONTINUE 
( 

C FINISH 
C 

C IY = DINT(DSEED) 

WRITE(13.50)IY 
50 FORMAT(UO) 

REWIND 13 

RETURN- 

END 
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c 

C SUBROUTINE NLCLAT 
C 

C THIS SUBROUTINE PRODUCES A CORRELATION MATRIX FROM NONLINEAR 
C TIME SEQUENCE IN AN ORDER WHICH IS COMPATIBLE WITH SUBROUTINE 
C SCHUR. 

C 

C WRITTEN 7 MAY 1985 
C 

q******************************** *************************************** 

C 

SUBROUTINE NLCLAT (Y,IYS,R,N) 

REAL*8 Y(10000),R(26,26),SUM,VEC(26) 

C 

C DEFINE CONSTANTS 
C 

MN = N*N 
MNPl = MN + 1 
IYSM2 = IYS - 2 
FIYSM2 = FLOAT(IYSM2) 

C 

C INITIALIZE R MATRIX TO ZERO 
C 

DO 20 I = l.MNPl 
DO 10 J = l.MNPl 
R(I.J) = 0.0 
10 CONTINUE 
20 CONTINUE 
C 

C BEGIN OUTER LOOP 
C 

DO 80 I = 3, IYS 

IR = 1 

VEC(IR) = Y(I) 

DO 50 MPl = l.N 
M0 - MPl - 1 
LLIM - 2* MPl - 1 
DO 40 L - 1 LLIM 
L0 L - 1 
11 Mo 
Jl = L0 2 

IF (MOD(L0.2).EQ.0) GO TO 30 

II Jl 
.11 = M0 

30 IR — IR — l 

YEC(IR) = COORD( Y(l-l).Il )*COOR D( Y(I-2),J 1 ) 

40 CONTINUE 

50 CONTINUE 

C 

C CALCULATE THE CORRELATIONS 
C 

DO 70 J l.MNPl 
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DO 60 K = J,MNPl 
R(J.K) = R(J.K) - VEC(J)*VEC(K) 

C WRITE(6.12)VEC(J),VEC(K).R(J.K) 

C12 FORMAT(3(2X,El2.5)) 

60 CONTINUE 

70 CONTINUE 
80 CONTINUE 
C 

C DIVIDE BY THE NUMBER OF DATA ELEMENTS CONSIDERED 
C 

DO 100 J = 1,MNP1 
DO 90 K = J,MNPl 

R(J,K) = R(J,K)/FIYSM2 
90 CONTINUE 
100 CONTINUE 
C 

C FILL IN THE SYMMETRIC HALF OF CORRELATION MATRIX 

C 

DO 120 I = 2,MNPl 
IMl = 1-1 
DO 110 J = l.IMl 
R(I.J) = R(J.I) 

110 CONTINUE 
120 CONTINUE 
C 
C 

RETURN 

END 
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q* ***** X ******************** X X * * X ± * * X * X X x * % * x. * * x % * X * * * * * * x * * x * * * * ± 

C 

C NONLINEAR WIENER MODELLING PROGRAM 
C THIS USES FUNCTIONS CONTAINED IN ROUTINE COORD 
C 

C MODIFIED 1 JAN 85, 14 JAN 85 
C 

q*********************************************************************** 

c 

DOUBLE PRECISION A(25,26),Z(25),X(15000),Y(15000),VAR,XA,YA,YHAT 
DOUBLE PRECISION ZZ(25),XMl,SEED,STORE(25) 

C 

C INPUT SIMULATION PARAMETERS 
C 

READ(13,77)IY 
REWIND 13 
WRITE(6,40) 

40 FORMAT(2X, ’MAGNITUDE OF NOISE’) 

READ(5,4l)VAR 

41 FORMAT(Fl2.5) 

WRITE(6,42) 

42 FORMAT(2X, ’NUMBER OF POINTS’) 

READ(5.43)N 

43 F0RMAT(I5) 

WRITE(6,44) 

44 FORM AT(2X. ’MAXIMUM POWER OF X ’) 

READ(5,45)LW 

45 FORMAT(U) 

WRITE(12.46)N 

WRITE(6.46)N 

46 FORM AT(2X. ’THE NUMBER OF POINTS WAS ’,15) 

WRITE) 12.4 7 )L\\ 

WRITE(6.47)LW 

47 FORMAT(2X.’THE MAXIMUM POWER OF X IS \Il) 
WRITE(6.48)YAR.VAR 

WRITE) 12.48) VAR. VAR 

48 FORMAT(2X.’THE NOISE IS UNIFORM FROM -’.F9.5.’ TO +’.F9.5) 

LW = LW - l 

VAR = YAR-2.0 
DO 132 1 = 1.25 
STORE(I) = 0.0 
132 CONTINUE 

X) 1 ) = YAR*(UR AND(IY) - .5) 

Y(l) 1.0 
DO 2 1 - 2 N 

X(I) YAR-(l'RAND)IY) - .5) 

Y(I) - L'NKNO W(X ( 1 ) .X) 1- 1 ) ) 

2 CONTINUE 

C A L L N C R L A T ( X . Y . A . LS Q W , L W . N ) 

C WRITE(6,134) 

C WRITE(12.134) 

Cl 34 FORMAT) 2X. ’NO INVERSION SOLUTION’) 

DO 142 I - 1. LSQW 
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ZZ(J) = A(J,LSQVV + 1) / A(J.J) 

STORE(J) = STORE(J) - ZZ(J)/25.0 
142 CONTINUE 
C DO 133 J = l.LW 
C IMIN = (J-l)'LVV i- 1 

C IMAX = J*LW 

C WRITE(6,173) (ZZ(I).I = IMIN, IMAX) 

C WRITE(12,173) (ZZ(I),I = IMIN, IMAX) 

C173 F0RMAT(5(2X,E12.5)) 

C133 CONTINUE 
C CALL SOLVE(A,Z,LSQW) 

C WRITE(6,135) 

C WRITE(12,135) 

C135 FORMAT(/2X, ’USING FULL MATRIX INVERSION’) 

C DO 27 J = l.LW 
C IMIN = (J-l)'LW + ] 

C IMAX = J*LW 

C WRITE(6,1 1) (Z(I),I = IMIN, IMAX) 

C WRITE) 12, 11) (Z(I),I = IMIN, IMAX) 

11 F0RMAT(5(2X,E12.5)) 

C27 CONTINUE 
131 CONTINUE 
WRITE(6.201 ) 

WRITE) 12,201) 

201 FORMAT) /2X. 'NO INVERSION SOLUTION AVERAGED OVER 25 RUNS’) 
DO 200 J = l.LW 
IMIN = (J - l)*L W + l 
IMAX = J * LW 

WRITE(6,11) (STORE(I).I = IMIN, IMAX) 

WRITE) 12.11) (STORE(I).I = IMIN.IMAX) 

200 CONTINUE 
WRITE(6,22) 

WR1TE(12.22) 

XMl = VAR*(URAND(1Y) - .5) 

DO 73 I = 1.10 

XA = VAR*(URAND(IY) - .5) 

YA UNKNO W (XA.XM 1 ) 

YHAT = 0.0 
YYHAT = 0.0 
DO 1 4 J — l.LW 
JMl .1-1 

X2 - COORD) XMl. JMl) 

DO 15 K l.LW 
KM 1 K - 1 
XI COORD(XA.KMl) 

C W RlTE(C.52)Xl.X2 

C 5 2 F O R M A T ( 2 ( 2 X . E 1 2 . 5 ) ) 

YHAT YHAT - Xl ’ X2 4 Z((J-1)*LW K) 

YYHAT YYHAT Xl 1 X2 * ZZ((J-1)’LW - K) 

15 CONTINUE 

14 CONTINUE 

ERROR (YA - YHAT) YA 
ZERROR (YA - YYHAT) YA 
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XMl = XA 

WRITE(6,17) YA,YHAT,YYHAT,ERROR.ZERROR 
17 FORMAT(5(2X.El2.5)) 

22 F0RMAT(/6X.’YA’, 12X,’YHATT0X. : YYHAT’,9X. ’ERROR’, 9X.’ZERROR’) 

WRITE(12,17) YA.YHAT,YYHAT.ERROR,ZERROR 
73 CONTINUE 
WRITE(13,77)IY 
77 FORMAT(IlO) 

STOP 

END 



oooooooo 



*****|!(;*?***x****i:i;****1.******»**x)i;**»*xix**x!kx»*i(;±*T**i****»**ii;****!([*x 



SUBROUTINE UNKNOW 

THIS ROUTINE DEFINES THE UNKNOWN SYSTEM 






DOUBLE PRECISION FUNCTION UNKNOW(X,Xl) 
DOUBLE PRECISION H(25),YHAT,XTl,XT2 
LW = 5 

H(l) = -2 

H(2) = -.4 
H(3) = .03 
H(4) = -.7 
H(5) = 0.0 
H(6) = .5 
H(7) = .35 
H(8) = .11 
H(9) = .9 
H(10) = 0.0 
H(ll) = .01 
H( 12) = 1.3 
H(13) = -.33 
H(14) = .7 
H( 1 5) = 0.0 
H( 1 6) = .43 
H(17) = .81 
H(18) = -.05 
H ( 1 9) = .4 
H(20) = 0.0 
H ( 2 1 ) = 0.0 
H (22 ) = 0.0 
H(23) = 0.0 
H(24) = 0.0 
H (25) = 0.0 
YHAT = 0.0 
DO 14 J - 1 .LW 
JM1 = .) - 1 

XT2 COORD(XlJMl) 

DO 15 K = l.LW 
KMl - K - 1 
XTl COORD(X.KMl) 

C URITE(G.52)XTl.XT2 

('52 FORMAT(2(2X El 2 5)) 

YHAT YHAT XT 1 ' XT2 ' H(JM1*LW - K) 
15 CONTINUE 
11 CONTINUE 

UN KNOW YHAT 

RETURN 

END 
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o n o o o o o o o o o o o o o o oooooooooo 



******************************** *x*****x************r**.*****x*x**xx-x*** 



FUNCTION COORD 

GENERATES OUTPUT OF THE FUNCTIONS BEING USED AS COORDINATES 
CREATED 23 AUG 84 

*********************************************************************** 



DOUBLE PRECISION FUNCTION COORD(X,I) 

USE LEGENDRE POLYNOMIALS 

IF (I.NE.O) GO TO 1 

Y = 1.0 
GO TO 30 

1 IF (I.NE.l) GO TO 2 

Y = 1.732051 *X 
GO TO 30 

2 IF (I.NE.2) GO TO 3 

Y - 3.354102*(X*X - l./S.) 

GO TO 30 

3 Y = 6.61438*(X*X*X - 3./5.*X) 

USE SIMPLE POWER SERIES TYPE POLINOMIALS 

Y = 1.0 

IF (I.EQ.0) GO TO 30 

Y = x **j 

30 COORD = DBLE(Y) 

RETURN 

END 



187 



o o 



c 

C SUBROUTINE SOLVE 

C 

C SUBROUTINE TO SOLVE A SYSTEM OF LINEAR EQUATIONS 
C USES GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING 

C 

q******************** ************************************ ************** 

c 

SUBROUTINE SOLVE(A,Z,K) 

DOUBLE PRECISION A(25, 26), Z(25),Y, TEMP, LARGE 

C 

KMl = K - 1 
KP1 = K + 1 



12 

C 



14 

c 

13 

C 



1 5 
11 
10 
c 
c 
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DO 10 L = l,KMl 
LP1 = L + 1 
DO 111= LPl.K 
LARGE = DABS(A(L,L)) 

LROW = L 
DO 12 M = LPl,K 

IF (DABS(A(M.Lj).LE.LARGE) GO TO 12 
LARGE - DABS(A(M,L)) 

LROW = M 
CONTINUE 

IF (L.EQ.LROW) GO TO 13 
DO 14 M = L,KPl 
TEMP - A(L.M) 

A(L.M) - A (LROW, M) 

A (LROW, M) - TEMP 
CONTINUE 

Y - A(I.L) 'A(L.L) 



DO 15 M = l.KPl 
A(I,M) = A ( I , M ) - A(L.M)*Y 
CONTINl E 
CONTINUE 
CONTINUE 



Z ( K ) - A(K.KPl) A(K.K) 
DO 16 L - I. KMl 
I K - L 
Zfl) A(I.KPl) 

KMl K - I 
DO 17 M l.K.MI 
J ~ K - M - 1 
Z(I) - Z(I) - A(1,J)*Z(J) 
CONTINUE 
m z(i) a ( i,i) 



188 



o o 



16 CONTINUE 



RETURN- 

END 



£X***slc*x*************x5fc*±****j|6x:t:****x******afcxx*x*xs|Ex*.xir*x*:fc:lc*slcjk****5|e:irsE 

C 

C NONLINEAR CORRELATION MATRIX CALCULATOR 
C 

C CREATED 23 AUG 84 

C USES FUNCTION UNKNOW TO DETERMINE FUNCTIONS FOR EXPANSION. 
C X - INPUT VECTOR 

C Y - OUTPUT VECTOR 

C PHI - CORRELATION MATRIX 

C - LAST COLUMN CONTAINS INPUT/OUTPUT CROSS-CORRELATION 
C 

C 

SUBROUTINE NCRLAT(X,Y,PHI.LSQW,LW,N) 

DOUBLE PRECISION PHI(25,26),X(1).Y(1),TX(5,5) 

C 

C 

NMl = N - 1 
LSQW — LW * LW 
LSQWPl = LSQW - 1 
C 
C 

DO 40 1=1. LSQW 
DO 39 J = l, LSQWPl 
PHI(LJ) = 0.0 

39 CONTINUE 

40 CONTINUE 
C 

C 

DO 2 M = 2.N 
MMl = M-l 
DO 3 I=1.LW 

IMl =1-1 

X2 = COORD(X(MMl ),IMl ) 

DO 4 J= l.LW 
JMl = .) - 1 

XI = COORD(X(.M),JM 1 ) 

C WR1TE(6,81 ) X1.X2 

('81 F0RMAT(2(2X.E12.5)) 

TX(l.J) = XI ’ X2 

C WRITE(6.80) X(K).X(MM 1 ).1.J,TX(I,J ) 

C80 FORM AT(2X,2(2X.El2.5).2(2X.I2).El 2.5) 

4 CONTINUE 

CONTINUE 

0 

C 

DO 5 1= l.LW 
DO 0 II ELW 
K (I-l)'LW II 

PHI(K. LSQWPl) PH1(K, LSQWPl) - Y(M)'TX(Lil) 

DO 7 J = 1.LW 
DO 8 JJ = 1,LW 

KK - (J-l)'LW - .1.1 
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PHI(K,KK) = PHl(K.KK) - TX(1,II)*TX(J,JJ) 
8 CONTINUE 

7 CONTINUE 

6 CONTINUE 

5 CONTINUE 

2 CONTINUE 
C 

DO 31 I=1,LSQW 
DO 32 J=l,LSQWPl 
PHI(I,J) = PHI(I,J)/FLOAT(NMl) 

32 CONTINUE 

31 CONTINUE 
C 
C 

C DO 17 I=1,LSQW 

C WRITE(11,16) (PHI(I.J), J = 1,LSQW) 

C16 FORMAT(/, 9(2X,E12.5)) 

C17 CONTINUE 

WRITE(11.19) (PHI(I,LSQWPl),I = l.LSQW) 

19 FORMAT(//, 9(2X,E12.5)) 

C 

C 

RETURN 

END 
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c 

C NONLINEAR WIENER MODELLING PROGRAM 

C THIS USES THE LMS ADAPTIVE ALGORITHM 14 FEB 84 (VALENTINE'S) 
C 

C UNKNOWN SYSTEM DEFINED IN FUNCTION UNKNOW 
C 

c 

DOUBLE PRECISION TX(5,5),H(5,5),VAR,X,XMl,Y,YHAT,SEED 

READ(13,77)IY 

REWIND 13 

WR1TE(6,40) 

40 FORMAT(2X, ’MAGNITUDE OF NOISE’) 

READ(5,41) VAR 

41 FORMAT(Fl2.5) 

WR1TE(6,42) 

4? FORMAT(2X, ’NUMBER OF POINTS’) 

READ(5,43)N 

43 FORMAT(I5) 

\VR1TE(6,44) 

44 FORMAT(2X. ’MAXIMUM POWER OF X ’) 

READ(5.45)LW 

45 FORMAT(Il) 

WRITE(6.7) 

7 FORMAT(2X, 'CONVERGENCE FACTOR’) 

READ(6.8)U 

8 F0RMAT(F12.5) 

WRITE(12,46)N 

WRITE(6.46)N 

46 FORMAT(2X,’THE NUMBER OF POINTS WAS ’.15) 

WRITE(12,47)LW 

WR1TE(6.47)LW 

47 FORMAT(2X,’THE MAXIMUM POWER OF X IS ’,11) 

W'R1TE(6.48) VAR.VAR 
WRITE(12.48)VAR.VAR 

48 FORM AT(2X. ’THE NOISE IS UNIFORM FROM -\F9.5,’ TO +’,F9.5) 
WRITE(6.49)U 

WRITEf 12.49)1’ 

49 FORMAT (2X.TIIE CONVERGENCE FACTOR WAS '.E12.5) 

LW LW ■ 1 

VAR = VAR‘2.0 

XMl - V \ R ’ (UR A N D(I Y) - .5) 

C 

( INITIALIZE THE II TENSOR 
(.' 

DO 10 I 1 LW 
DO 9 .) 1 LW 

ll(I.J) - 0.0 

9 CONTINUE 

10 CONTINUE 
C 

C BEGIN THE ITERATION LOOP 
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c 

DO 2 K = 2,N 

X = VAR* (UR AND (IY) - .5) 

Y = UNKNO W(X.XMl) 

C WRITE(6.97)I,X,Y 

C97 F0RMAT(2X.I5.2X,E12.5,2X.E12.5) 

YHAT = 0.0 
C 

C CALCULATE THE MEASUREMENT TENSOR AND THE OUPUT ESTIMATE 
C 

DO 14 I = 1,LW 
IMl = 1-1 

X2 = COORD(XMl,IMl) 

DO 15 J = 1,LW 
JMl = J - 1 
XI = COORD(X,JMl) 

C WRITE(6,52)Xl,X2 

C52 FORM AT(2(2X.El 2.5)) 

TX(I.J) = Xl*X2 

YHAT = YHAT + TX(I,J) * H(I,J) 

15 CONTINUE 

14 CONTINUE 
C 

C CALCULATE THE NEW VALUE OF THE TENSOR 
C 

ERROR = Y - YHAT 
DO 5 I = l.LW 
DO 4 J = l.LW 

H(I.J) = H(I,J) + 2.0*U*ERROR*TX(I,J) 

4 CONTINUE 

5 CONTINUE 
XMl = X 

C 

C PRINT THE NEW VALUE OF THE TENSOR 
C 

KMl = K - 1 

IF (KMl. NE. (KMl 100)*100)CO TO 2 

WRITE(6.13)KM1 

WRITE) 1 2.1 3)KMl 

13 FORMAT) 2X. ITERATION NUMBER .15) 

WR1TE(C 135) 

W RITE) 12.135) 

135 FORMAT(2X 'THE NEW VALUE OF THE H TENSOR IS ) 

DO 27 I l.LW 

\YRITE(6.1I) (H(L.I).J l.LW ) 

W RITE( 12.11) (H(I.J).J l.LW ) 

II FORM AT(5(2X.El 2.5)) 

27 CONTINUE 

W RITE(6.22) 

WRITEf 12.22) 

WRITE(6.17) Y. YHAT. ERROR 
17 FORM AT(5(2X.El 2.5)) 

22 FORMAT(6X. ! Y . 12X,' YHAT’ 1 OX. ERROR') 
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WRITE(12,17) Y.YHAT, ERROR 
CONTINUE 
\VRITE(13,77)IY 
FORMAT(IlO) 

STOP 

END 
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